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Abstract 

Estuarine  density  stratification  may  be  controlled  primarily  by  cross-shore  processes 
(analogous  to  longitudinal  control  in  narrow  estuaries),  or  by  both  cross-  and  alongshore 
processes  (typical  of  coastal  plumes).  Here  field  observations  and  numerical  modeling  are 
used  to  investigate  stratification  on  the  low-sloped,  periodically  inundated  Skagit  Bay 
tidal  flats.  Advection  of  stratification  by  the  depth-averaged  velocity,  straining  of  the 
horizontal  density  gradient  by  velocity  shear,  and  turbulent  mixing  are  shown  to  be  the 
dominant  processes.  On  the  south-central  flats  (near  the  south  fork  river  mouth)  velocities 
are  roughly  rectilinear,  and  the  largest  terms  are  in  the  major  velocity  direction  (roughly 
cross-shore).  However,  on  the  north  flats  (near  the  north  fork  river  mouth),  velocity 
ellipses  are  nearly  circular  owing  to  strong  alongshore  tidal  flows  and  alongshore 
stratification  processes  are  important.  Stratification  was  largest  in  areas  where  velocities 
and  density  gradients  were  aligned.  The  maximum  stratification  occurred  during  the 
prolonged  high  water  of  nearly  diurnal  tides  when  advection  and  straining  with  relatively 
weak  flows  increased  stratification  with  little  mixing.  Simulations  suggest  that  the 
dominance  of  straining  (increasing  stratification)  or  mixing  (decreasing  stratification)  on 
ebb  tides  depends  on  the  instantaneous  Simpson  number  being  above  or  below  unity. 

Thesis  Supervisor:  Britt  Raubenheimer 

Title:  Associate  Scientist  in  Applied  Ocean  Physics  and  Engineering 
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1  Introduction 


1.1  Overview 

Shallow,  periodically  submerged,  low-sloped  tidal  flats  bordering  coasts  and  estuaries 
may  affect  water  circulation  and  temperature  in  the  surrounding  region.  They  can  be 
protective  barriers  to  low-lying  land  as  well  as  sources  of  economic  activity,  habitats  for 
juvenile  fish,  and  systems  of  natural  water  purification  by  the  benthic  organisms  living 
inside  the  flats.  The  viability  of  tidal  flat  systems  near  river  mouths  is  affected  by  the 
density  distribution,  including  stratification,  and  by  circulation  and  sediment  transport, 
which  also  are  affected  by  stratification.  However,  most  prior  studies  of  density 
stratification  near  river  mouths  have  focused  on  estuaries  and  coastal  plumes,  and  there 
have  been  few  studies  of  stratification  on  deltaic  tidal  flats.  It  is  expected  that  mixing  and 
stratifying  processes  may  be  different  on  flats  owing  to  the  shallow  water  depths.  Here, 
observations  and  model  predictions  will  be  used  to  examine  the  processes  controlling  the 
generation  and  destruction  of  stratification  on  tidal  flats,  leading  to  a  better  understanding 
of  these  systems,  which  may  allow  improved  predictions  of  the  effects  of  human  activity. 

The  specific  goals  of  this  work  are  to: 

•  Determine  the  dominant  processes  controlling  changes  in  stratification  over  a  tidal 
cycle. 

•  Investigate  the  temporal  changes  in  the  stratifying  processes  using  field  data. 

•  Use  model  output  to  examine  the  spatial  variability  of  the  different  processes. 

•  Examine  the  relative  importance  of  cross-  and  alongshore  processes. 

•  Relate  spatial  and  temporal  differences  in  stratification  to  tidal  and  estuarine 
parameters. 

•  Indicate  effects  of  stratification  on  circulation. 
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1.2  Background 


Tidal  flats  with  broad,  low-slope  intertidal  (intermittently  submerged)  shoals,  or  “flats”, 
traversed  or  bordered  by  deeper  riverine  or  tidal  channels  are  common  along  the  edges  of 
estuaries.  The  channels  are  significant,  being  the  primary  conduits  of  freshwater  and 
sediment  to  the  coast  (Chen  et  al.,  2010),  but  the  flats  contain  the  majority  of  area  and 
habitat. 

Although  tidal  flats  form  a  relatively  small  portion  of  the  total  coastal  and  estuarine  area, 
they  can  affect  the  circulation  patterns  throughout  larger  basins  by  increasing  heating 
(Kim  et  al.,  2010;  Kim  and  Cho,  201 1),  bed  friction  (Nicolle  and  Karpytchev,  2007),  and 
tidal  volume  storage  (Friedrichs  and  Aubrey,  1988),  and  by  altering  basin  resonance 
characteristics  (Fortunato  et  al.,  1997,  1999)  and  nonlinear  tidal  interactions  (Speer  and 
Aubrey,  1985;  Fortunato  et  al.,  1999;  Blanton  et  al.,  2002).  Tidal  flats  also  can  provide  a 
coastal  buffer  (Kirby,  2000;  Kim,  2003)  and  important  habitat  for  fish  and  game 
(Grossman  et  al.,  2007).  A  healthy  wetland  ecosystem  will  contain  intertidal  flats  as  well 
as  marshes  and  subtidal  (always  submerged)  areas  (Havens  et  al.,  1995).  The  benthic 
organisms  in  the  flats  constitute  a  natural  system  of  water  purification,  removing 
suspended  organic  matter  that  can  be  detrimental  to  local  fisheries  (Suzuki,  2001).  Some 
of  the  work  relating  to  tidal  flats  has  focused  on  the  hydrodynamics  in  the  deeper 
channels  (e.g.,  Ralston  and  Stacey  2005a, b,  2007),  sediment  transport  (e.g.  Chen  et  al., 
2010;  Talke  and  Stacey  2008;  Lee  et  al.,  2004),  or  the  effects  of  the  flats  on  their 
surroundings  (e.g.  Kim  et  al.,  2010;  Nicolle  and  Karpytchev,  2007;  Fortunato  et  al., 

1999).  However,  understanding  the  complicated  circulation  and  salinity  characteristics 
resulting  from  the  shallow  depths,  periodic  inundation,  tidal  flows,  and  river  discharge  on 
the  shoals  is  vital  to  managing  these  natural  resources. 

Stratification,  the  layering  of  different  densities  of  water,  affects  the  location  and 
magnitude  of  fluid  turbulence,  which  affects  the  flow  field  and  the  sediment  transport.  In 
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particular,  stratification  affects  the  viscous  term  in  the  momentum  balance,  which 
changes  the  velocity  profile  by  changing  the  production  and  magnitude  of  turbulence.  For 
example,  a  well-mixed  water  column  over  a  rough  bed  has  a  shear  flow  with  a  log  profile 
given  by  —  log(z/z0),  where  zo  is  the  roughness  length,  z  is  height  above  bed,  u*  is  the 

K 


friction  velocity,  and  k  is  von  Karman’s  constant  (von  Karman,  1930).  For  the  same 
external  forcing,  but  with  stable  stratification,  the  flow  profile  will  be  log-linear  and  goes 

I#  /  z  \ 

as  —  log  I — I-  bz ),  where  b  is  a  scaling  parameter  that  depends  on  the  stratification. 

rC  \Zq  / 


Under  stably  stratified  conditions,  the  flows  that  are  far  from  the  bed  are  larger  than  they 
would  be  given  a  log  profile  because  some  of  the  turbulent  energy  is  dissipated  by 
destabilizing  the  fluid  rather  than  by  slowing  the  upper  water  column  (Turner,  1973). 
Overall,  stratification  tends  to  decouple  the  water  near  the  bed  from  that  near  the  surface. 
However,  although  stratification  tends  to  suppress  turbulent  production  near  the  bed, 
velocity  shear  can  produce  instabilities  at  the  pycnocline  when  the  Richardson  number 
(ratio  of  buoyancy  to  shear)  is  small  (Ralston  et  ah,  2010b;  Peters  and  Bockhorst,  2000; 
Peters  1997). 


Stratification  can  affect  sediment  transport  through  the  effects  on  the  flow  profile  and 
turbulence.  In  a  stratified  fluid,  reduced  turbulence  leads  to  reduced  re-suspension,  and 
thus  the  sediment  is  more  likely  to  settle  out  and  collect  on  the  bed  (Dyer,  1986).  This 
phenomenon  enhances  the  concentration  of  sediment  in  an  estuarine  turbidity  maximum 
(Geyer,  1993).  The  effect  is  most  pronounced  on  silt-sized  particles.  Substantially  coarser 
particles  fall  out  of  suspension  too  easily  for  the  stratification  to  have  much  of  an  effect, 
and  very  fine  particles  settle  slowly  enough  that  they  remain  in  suspension  with  or 
without  stratification  or  additional  turbulence. 


Intra-  and  intertidal  (within  and  averaged  over  the  tidal  cycle,  respectively)  temporal 
variations  of  stratification  in  estuaries  and  on  estuarine  tidal  flats  have  been  shown  to 
change  circulation  patterns  (Monismith  et  al.,  1996;  Geyer  et  al.,  2000;  Stacey  et  al., 
2001;  Ralston  and  Stacey,  2005a;  Becker  et  al.,  2009;  Cheng  et  al.,  2009),  suppress 


9 


turbulence  (Nepf  and  Geyer,  1996;  Peters  and  Bokhorst,  2001;  Rippeth  et  al.,  2001; 
Ralston  and  Stacey,  2005b;  Stacey  et  al.,  1999;  Ralston  et  al.,  2010b;  Wang  et  al.,  2011), 
and  reduce  bottom  stress  and  suspended  sediment  concentration  (Chant  and  Stoner,  2001; 
Ralston  and  Stacey,  2007;  Chen  et  al.,  2010). 

The  potential  energy  anomaly  (Simpson  and  Bowers,  1981),  the  amount  of  energy  per 
unit  volume  required  to  homogenize  the  water  column,  often  is  used  to  quantify  changes 
in  stratification  (Simpson  et  al.,  1990;  Wiles  et  al.,  2006;  Marques  et  al.,  2010,  2011; 
Ralston  et  al.,  2010b).  The  potential  energy  anomaly,  denoted  by  d>,  is  defined  as: 

$  =  ~p)zdz  (1) 

where  g=9.8  m/s  is  gravitational  acceleration,  D=ri+h  is  the  water  depth,  r\  is  the  surface 
elevation,  h  is  bed  level  below  the  mean  surface,  p  is  the  water  density,  the  overbar 
denotes  a  depth-average,  and  z  is  the  vertical  coordinate,  positive  up. 

In  many  estuaries  the  potential  energy  anomaly  owing  to  longitudinal  tidal  straining  (the 
vertically  sheared  velocity  profile  acting  on  the  along-channel  horizontal  density 
gradient)  increases  stratification  and  suppresses  turbulent  mixing  on  the  ebb  and 
decreases  stratification  and  enhances  mixing  on  the  flood  (Simpson  et  al.,  1990;  Chant 
and  Stoner,  2001;  Burchard  and  Hofmeister,  2008).  In  shallow  salt- wedge  estuaries 
longitudinal  advection  may  enhance  straining  effects  on  the  flood  and  may  oppose  the 
straining-induced  increase  in  stratification  on  the  ebb  (Giddings  et  al.  2011).  Vertical 
advection  processes  also  may  be  important  in  regions  with  large,  spatially 
inhomogeneous  horizontal  density  gradients  (Nepf  and  Geyer,  1996;  Burchard  and 
Hofmeister,  2008;  de  Boer  et  al.,  2008;  Marques  et  al.,  2011),  whereas  in  estuaries  with 
complex  bathymetry  and  in  coastal  regions  near  river  plumes  (e.g.,  regions  of  freshwater 
influence  or  ROFIs)  both  cross-  and  alongshore  processes  contribute  to  the  stratification 
(Valle-Levinson  and  Atkinson,  1999;  Lacy  et  al.,  2003;  de  Boer  et  al.,  2008;  Marques  et 
al.,  2010).  Furthermore,  nonlinear  effects  may  be  important  near  river  mouths  and 
estuarine  inlets  (de  Boer  et  al.,  2008;  Marques  et  al.,  2010),  and  wind-driven  currents  can 
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affect  the  potential  energy  anomaly  balance  in  estuaries  and  ROFIs  during  storms  (Yang 
and  Khangaonkar,  2009;  Marques  et  al.,  2010,  2011). 

Overall,  the  magnitudes  of  the  river  discharge,  tidal  and  wind-driven  currents,  and 
horizontal  density  gradients  influence  which  processes  dominate  the  stratification  balance 
(Nepf  and  Geyer,  1996;  Burchard  and  Hofmeister,  2008;  de  Boer  et  ah,  2008;  Hofmeister 
et  ah,  2009).  These  forcing  mechanisms  and  the  stratification  can  vary  on  seasonal 
(Marques  et  ah,  2010),  spring-neap  (Peters,  1997),  storm  (Marques  et  ah,  2010,  2011), 
and  tidal  timescales  (Simpson  et  ah,  1990;  Nepf  and  Geyer,  1996;  Stacey  et  ah,  1999; 
Rippeth  et  ah,  2001).  In  the  Merrimack  (Massachusetts)  during  high  river  discharge 
maximum  stratification  occurred  on  the  flood  and  the  dominant  processes  were  the 
advection  and  mixing  of  a  highly  stratified  salt  wedge  (Ralston  et  ah,  2010a).  During 
more  moderate  river  discharge,  maximum  stratification  occurred  on  the  ebb  and  tidal 
straining  and  mixing  were  the  dominant  processes.  In  the  Snohomish  estuary 
(Washington)  the  dominant  processes  were  the  upstream  advection  of  a  salt  wedge, 
straining  of  the  stratified  water  by  the  sheared  flows,  and  mixing  (Giddings  et  ah,  201 1). 
During  spring  tides,  strong  flood  flows  mixed  the  water  column  offshore  of  (behind)  the 
salt  wedge  front,  and  thus  stratification  first  increased  as  the  salt  wedge  was  advected 
upstream  and  then  decreased  as  the  front  passed.  In  contrast,  during  neap  tides  flood 
flows  (and  mixing)  were  weaker,  and  thus  the  water  behind  the  salt  wedge  remained 
stratified. 

The  stratification  balance  and  the  dominant  processes  are  spatially  variable.  For  example, 
numerical  simulations  of  a  region-of- freshwater-influence  (ROFI,  a  coastal  river  plume) 
suggested  that  the  dominant  processes  change  with  distance  from  the  freshwater  source. 
Alongshore  advection  was  found  to  dominate  at  the  edge  of  the  bulge  but  cross-shore 
straining  dominated  farther  down  stream  in  the  coastal  current  region  (de  Boer  et  al., 
2008)  despite  strong  alongshore  velocities.  Simulations  of  a  salt-wedge  channel  estuary 
showed  a  balance  of  advection,  straining,  and  vertical  processes  at  the  salt  wedge  front 
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but  strain-induced  periodic  stratification  farther  downstream  (Burchard  and  Hofmeister, 
2008).  Just  as  stratification  can  be  different  between  in  the  thalweg  (main  channel)  and 
over  the  shoals  of  estuaries  (Cheng  et  al.,  2009)  because  there  are  additional  layers  of 
dense  water  at  the  bottom  of  the  deeper  area,  stratification  on  the  shoals  of  estuarine  tidal- 
flat  regions,  which  make  up  the  majority  of  the  total  surface  area,  can  differ  significantly 
from  that  in  the  river  and  tidal  channels.  In  some  cases  the  water  over  the  shoals  may 
remain  well  mixed  and  saline  despite  stratifying  fresh-water  influence  in  the  channels 
(Ralston  and  Stacey,  2005b).  Stratification  can  result  in  or  near  channels  from  density 
gradients  between  channels  and  shoals  (Ralston  and  Stacey,  2005a)  by  having  fast  flows 
bring  different  density  water  down  the  channel  than  is  present  over  the  shoals. 

Similar  to  macrotidal  salt-wedge  estuaries  on  the  Merrimack  (Massachusetts)  (Ralston  et 
al.,  20 10a, b),  Columbia  ( Washington/ Oregon  border)  (Jay  and  Smith,  1990),  Fraser 
(British  Columbia)  (Geyer  and  Farmer,  1989),  and  Snohomish  (Washington)  rivers 
(Wang  et  al.,  2009;  Giddings  et  al.,  2011)  where  the  length  of  the  salinity  transition  is 
similar  to  the  tidal  intrusion,  tidal  flats  and  the  associated  channels  have  strong  cross¬ 
shore  density  gradients,  and  periodically  are  strongly  stratified.  Tidal-flat  channels  can  be 
similar  to  narrow  estuaries,  especially  at  low  tide  and  during  periods  of  high  river  runoff. 
For  example,  stratification  increasing  on  ebbs  and  decreasing  on  floods  in  a  San 
Francisco  Bay  (California)  tidal-flat  channel  during  the  spring  freshet  was  primarily  a 
result  of  straining  of  the  longitudinal  density  gradient  and  advection  of  the  salinity  front 
(Ralston  and  Stacey,  2005a, b,  2007).  However  on  tidal  flats  the  tidal  range  is  greater  than 
the  mean  water  depth.  Thus  the  effects  of  depth  changes  can  be  significant  (Ralston  and 
Stacey,  2005b,  2007;  Giddings  et  al.,  2011),  and  the  effects  of  the  periodic  inundation 
and  drying  of  tidal  flats  on  the  balance  of  stratification  is  uncertain. 

Water  levels,  velocities,  and  densities  in  estuaries  (Ralston  et  al.,  20 10a, b)  and  on  tidal 
flats  (Yang  and  Khangaonkar,  2009;  Kim  and  Cho,  2011;  Ralston  et  al.,  2012),  have  been 
predicted  with  the  Finite  Volume  Coastal  Ocean  Model  (FVCOM),  a  three-dimensional 
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finite  volume,  unstructured  grid,  primitive  equation  model  (Chen  et  al.,  2003)  that  was 
developed  to  simulate  coastal  regions.  Model  comparisons  have  shown  that  tidal 
amplitudes  and  phases  in  regions  with  irregular  coastlines  or  bathymetry  are  predicted 
significantly  better  by  models  with  unstructured  grids  than  by  models  with  rectangular 
mesh  elements  (Chen  et  al.,  2003).  In  addition,  the  wet/dry  point  technique  used  in 
FVCOM  enables  accurate  predictions  in  areas  that  are  periodically  submerged  or 
exposed.  For  the  Baeksu  tidal  flats  (South  Korea),  tidal  elevation  predictions  from 
FVCOM  agreed  well  with  observations  (r  >0.95)  while  tidal  velocity  ellipse  predictions 
were  of  similar  magnitude  but  slightly  rotated  counterclockwise  from  the  observations 
(Kim  and  Cho,  2011).  Observations  of  water  level,  velocity,  and  salinity  have  been 
predicted  well  (correlation  r  >  0.9)  in  estuaries  (e.g.,  Ralston  et  al.,  2010a)  and  on  tidal 
flats  (Kim  and  Cho,  2011;  Ralston  et  al.,  2012).  Tidal  evolution  of  salinity  fronts  also  is 
reproduced,  although  time  series  of  surface  salinity  are  more  difficult  to  simulate  owing 
to  large  spatial  and  temporal  variations  (up  to  20  PSU  over  100  m)  near  fronts. 

In  this  thesis,  observations  and  FVCOM-simulations  of  water  levels,  velocities,  and 
densities  in  Skagit  Bay,  Washington  (described  further  in  Chapter  2)  are  used  to 
investigate  the  processes  leading  to  production,  destruction,  and  movement  of 
stratification  on  tidal  flats.  Overall,  the  following  questions  are  addressed  in  Chapters  3 
and  4: 

1 .  What  processes  affect  the  stratification,  and  what  causes  spatial  variations  in  the 
magnitudes  of  these  processes? 

2.  What  is  the  relative  importance  of  cross-  and  alongshore  processes,  and  what 
controls  this  balance? 

3.  Where  is  stratification  generated  and  destroyed? 

In  Chapter  3,  the  observations  collected  during  August  2009  on  the  north  flats  are  used  to 
examine  the  relative  importance  of  across-  and  alongshore  straining,  advection,  and 
mixing  in  a  region  with  strong  alongshore  flows.  Prior  studies  have  shown  that  FVCOM 
accurately  predicts  water  depths  and  velocities  near  the  deepest  river  channels,  and 
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salinity  in  the  tidal  channel  along  the  western  edge  of  the  bay  (Yang  and  Khangaonkar, 
2009).  The  model  also  has  been  shown  to  reproduce  the  water  levels,  velocities,  and 
salinities  in  a  tidal  flat  channel  near  the  south  fork  of  the  Skagit  River  (Ralston  et  al., 
2012).  In  Chapter  4,  this  field- verified  model  is  used  to  examine  the  spatial  dependence 
of  the  processes  affecting  stratification.  The  flats  extend  from  the  mouth  of  the  north  fork 
of  the  Skagit  River,  about  10  km  south  east  to  the  mouth  of  the  south  fork  of  the  Skagit 
River.  Freshwater  discharge  is  large  near  the  river  mouths  at  either  end  of  the  bay,  but 
there  is  little  direct  freshwater  discharge  from  the  marshes  between  the  two  river  mouths. 
Thus,  the  processes  affecting  stratification  on  the  flats  are  examined  as  a  function  of 
proximity  to  the  freshwater  supplies.  Furthermore,  the  tidal  velocities  are  nearly 
rectilinear  and  aligned  with  the  density  gradients  on  the  south-central  flats  (similar  to  an 
estuarine  channel),  but  have  similar  cross-  and  alongshore  magnitudes  on  the  north  flats 
(similar  to  a  ROFI),  so  the  stratification  processes  for  these  different  types  of  coastal 
regimes  are  compared  and  contrasted. 
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2  Site  Description 


2.1  Geographic  setting 

The  Skagit  Bay  tidal  flats,  near  La  Conner,  WA,  have  an  area  of  about  100  km  (Fig.  2- 
1).  A  deep  channel  runs  along  the  edge  of  Whidbey  Island,  which  forms  the  western 
border  of  Skagit  Bay.  To  the  north  and  south,  Skagit  Bay  connects  with  the  Strait  of  Juan 
de  Fuca  and  the  rest  of  Puget  Sound  via  Deception  Pass  and  Saratoga  Passage, 
respectively.  Tides  propagate  northwestward  from  Saratoga  Passage  towards  Deception 
Pass. 

About  5  km  upstream  of  the  flats,  the  Skagit  River  splits  into  north  and  south  forks  that 
carry  approximately  60%  and  40%  of  the  flow,  respectively  (Grossman  et  al.,  2007;  Yang 
and  Khangaonkar,  2009).  Data  from  the  field  study  discussed  in  Chapter  3  were  collected 
on  the  shoals  between  the  north  and  south  forks  during  the  month  of  August,  2009. 
Additional  data  collected  in  July  are  compared  with  a  numerical  model  in  Chapter  4. 
Typical  monthly-average  river  discharge,  as  measured  at  Mt.  Vernon  (upstream  of  the 
fork)  is  about  500  m  /s  with  a  maximum  in  June  during  the  spring  freshet  averaging  700 
m  /s  and  a  minimum  at  the  end  of  the  summer  in  September  averaging  about  250  m  /s. 
The  field  study  period  of  early  July  to  late  August  was  selected  to  capture  as  much 
variability  in  river  discharge  as  possible  over  a  short  field  season.  River  flow  during  the 
summer  of  2009  was  lower  than  average,  about  200  m  /s  during  the  August  study  period 
and  between  300  and  500  m3/s  during  the  July  period  (USGS  gage  12200500 
http://waterdata.usgs.gov/nwis/nwisman/7site  no=12200500).  Numerous  small  channels 
(depth  0( 0.10-0.25  mj)  split  off  from  the  north  fork  of  the  Skagit  and  extend  across  the 
marshes  onto  the  tidal  flats  near  the  measurement  locations  (Elgar  and  Raubenheimer, 
2011;  Webster  et  al.,  2012).  However,  the  majority  of  the  discharge  from  the  north 
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(south)  fork  exits  onto  the  flats  about  2.5  km  northwest  (4.5  km  southeast)  of  the  field 
study  area. 

In  contrast  with  many  intertidal  flats  in  low  energy  environments  (Banas  et  al.,  2004;  Lee 
et  al.,  2004;  Fan  et  al.,  2006),  the  sediment  on  the  Skagit  flats  is  primarily  sandy  (Webster 
et  al.,  2012).  The  cross-shore  bed  slope  on  the  flats  ranges  between  approximately  1/2000 
and  1/1000.  The  spring  tidal  range  is  about  4  m,  and  much  of  the  lower  part  of  the  flats 
are  dry  at  lower  low  tide.  Thus,  the  tidal  range  is  larger  than  the  mean  water  depth,  in 
contrast  to  deeper  estuaries  and  ROFIs.  Tides  are  mixed  (Fig.  2-2),  and  can  be  nearly 
semidiurnal  (denoted  here  as  “type  1”)  or  nearly  diurnal  (“type  2”). 

The  two  tide  types  were  defined  according  to  the  diurnal  inequality  to  enable  averaging 
over  tidal  cycles.  For  the  field  data  discussed  in  Chapter  3,  cycles  where  the  difference 
between  the  heights  of  the  two  low  tides  measured  on  the  mid  flats  (where  the  bed 
sometimes  is  dry)  in  one  diurnal  period  was  <  1/4  of  the  local  daily  tidal  range  were 
defined  as  type  1  tides,  and  cycles  where  the  difference  was  >  1/3  of  the  tidal  range  were 
defined  as  type  2.  For  the  model  results  (Chapter  4),  the  criteria  for  the  difference  in 
offshore  (always  submerged)  low  tide  levels  were  <5/12  and  >  7/12  of  the  daily  range 
for  type  1  and  type  2  tides,  respectively.  The  use  of  offshore  rather  than  local  tide  levels 
necessitated  the  change  in  empirical  criteria.  The  fortnightly  transition  between  low  and 
high  diurnal  inequality  with  the  latter  manifesting  as  a  long  high  tide  during  type  2  tides 
is  common  along  the  west  coast  of  North  America,  and  in  particular  the  local  region 
(Nidzieko  2010). 
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Figure  2-1:  (a)  Map  of  Skagit  Bay  and  surrounding  area  (from  NOAA/NOS  Medium  Resolution 
Coastline  Database)  and  instrument  array  (black  circles  and  yellow  star  indicate  where  density 
profiles  and  co-located  velocity  and  density  profiles,  respectively,  were  obtained).  Bottom 
pressure  (and  water  depth)  was  measured  at  all  instrument  locations.  Shading  indicates  four  levels 
of  bathymetry  with  darker  being  deeper  water.  The  smaller  map  (NOAA  World  Vector  Shoreline 
Database)  on  the  right  shows  the  location  of  Skagit  Bay  on  the  Pacific  Northwest  coast.  Positive 
cross-shore  is  toward  the  northeast  and  positive  alongshore  is  toward  the  northwest,  (b) 
Bathymetry  (relative  to  NAVD88)  around  the  instrument  area  rotated  into  alongshore  and  cross¬ 
shore  directions. 
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2.2  Modeled  circulation  and  salinity  distributions 


Differences  in  mean  water  level  and  tidal  range  between  Deception  Pass  and  Saratoga 
Passage  result  in  strong  flows  through  Skagit  Bay,  particularly  in  the  deep  channel 
(“gutter”)  on  the  western  edge  of  the  bay.  The  circulation  and  surface  salinity  fields 
during  different  tidal  stages  are  illustrated  using  snapshots  of  the  FVCOM  numerical 
simulations.  Northing  and  Easting  are  for  UTM  zone  10U. 


Time  (day  in  August  2009) 


Figure  2-2:  (a)  Recorded  water  levels  in  Seattle,  WA  relative  to  MLLW  (from  NOAA  tidal  gage 
#9447130)  (b)  Measured  water  depth  versus  time  at  the  central  sensor  location  with  type  1  and 
type  2  tides  shaded  in  yellow  and  purple,  respectively 


During  the  type  1  strong  flood,  water  spreads  out  of  the  gutter  and  flows  onto  the  flats 
towards  the  east  on  the  south  flats  and  towards  the  northeast  on  the  north  flats  (Fig.  2-3a). 
At  this  mid-flood  stage,  the  fresh  surface  water  (remaining  from  the  last  ebb)  has  been 
transported  shoreward  towards  the  marshes  with  more  saline  surface  water  behind.  At 
high  tide,  flow  magnitudes  are  similar  to  those  on  flood  and  ebb  and  are  directed  towards 
Deception  Pass  except  on  the  south-central  flats,  which  are  shadowed  by  Camano  Island, 
which  forms  the  southern  border  of  the  tidal  flats  (Fig.  2-3c).  Saline  water  has  propagated 
onshore  and  a  roughly  alongshore -homogeneous  sharp  surface  front  occurs  on  the  upper 
flats  (Fig.  2-3d). 
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During  the  type  1  ebb,  the  water  flows  towards  the  west  over  most  of  the  flats,  and 
towards  Deception  Pass  in  the  northern  part  of  the  gutter  (Fig.  2-3e).  The  sharp  surface 
salinity  front  is  transported  offshore,  with  north  fork  discharge  rounding  the  islands  and 
going  northwest  towards  Deception  Pass  (Fig.  2-3f).  Near  low  tide,  when  water  is 
shallow  on  the  flats,  flows  on  the  south  and  central  flats  are  small,  flows  on  the  north  flats 
are  slightly  larger  and  southeastward,  and  flows  in  the  gutter  are  strong  with  water 
flowing  toward  Saratoga  Passage  (Fig.  2-3g).  Strong  flows  also  occur  in  the  south  fork 
channels  as  water  drains  from  the  river.  The  surface  water  remaining  on  the  flats  and  in 
the  gutter  is  brackish  or  fresh.  Throughout  the  tide  cycle,  the  freshest  surface  water  is 
contained  near  the  south  fork  channels  and  near  the  mouth  of  the  north  fork,  with 
partially  mixed  water  along  the  marshes  between  the  river  mouths  (Fig  2-3b,  d,  f,  and  h, 
regions  of  dark  blue). 

Flows  and  surface  salinities  during  the  type  2  strong  flood,  the  following  (first)  high  tide, 
and  the  strong  ebb  (not  shown)  are  similar  to  those  during  type  1  tides.  During  the  type  2 
weak  ebb  (not  shown),  flows  are  weaker  than,  but  similar  to,  those  during  the  type  1 
strong  ebb,  and  the  surface  freshwater  propagates  offshore  slowly.  However,  in  contrast 
to  the  type  1  low  tide,  during  the  type  2  higher  low  water  the  majority  of  the  flats  remain 
submerged  and  the  southeastward  flows  occur  over  most  of  the  north  flats  (similar  to  the 
type  1  low  tide,  flows  are  small  on  the  south  flats)  (compare  Fig.  2-4a  with  Fig.  2-3g). 
Furthermore,  strong  alongshore-inhomogeneous  surface  salinity  fronts  occur  near  the 
offshore  edge  of  the  flats  and  over  the  north-central  flats.  During  the  type  2  weak  flood 
flows  are  small  (Fig.  2-4c)  and  form  a  counterclockwise  circulation  cell  just  south  of  the 
north  fork,  which  results  in  spreading  of  the  freshwater  plume  (Fig.  2-4d).  Flows  and 
surface  salinities  during  the  second  type  2  high  tide  (Fig.  2-4e  and  f)  are  similar  to  those 
during  the  type  1  high  tide  (Fig.  2-3c  and  d),  but  the  fresh  surface  water  extends  farther 
from  the  marshes. 
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Figure  2-3:  Snapshots  of  depth-averaged  velocity  (a,  c,  e,  and  g)  and  surface  salinity  (b,  d,  f,  and 
h)  for  type  1  flood  (a  and  b),  high  water  (c  and  d),  ebb  (e  and  f),  and  low  water  (g  and  h) 
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Figure  2-4  Snapshots  of  depth-averaged  velocity  (a,  c,  and  e)  and  surface  salinity  (b,  d,  and  f)  for 
type  2  higher  low  water  (a  and  b),  weak  flood  (c  and  d),  and  high  water  after  the  weak  flood  (e 
and  f). 


The  tides  have  been  separated  into  type  1  and  type  2  to  account  for  the  temporally 
dependent  structure  observed  in  this  area.  Type  1  tides  are  similar  to  the  semi-diurnal 
tides  often  observed  along  the  U.S.  east  coast.  In  contrast,  the  type  2  tides  have  an 
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extended  period  of  high  water  during  the  weak  ebb  and  flood.  The  weak  ebb  and  flood 
have  distinct  circulation  patterns,  and  consequently  the  salinity  transport  and  the 
evolution  of  stratification  can  be  expected  to  differ  from  the  stronger  floods  and  ebbs. 

The  weak  flows  during  the  weak  flood  and  ebb  also  allow  density  driven  flows  to  become 
more  prominent,  which  is  shown  to  cause  differences  in  stratification  between  the  two 
types  of  tides. 
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3  Observations 


This  chapter  is  based  on  a  paper  submitted  to  Continental  Shelf  Research  entitled 
Processes  controlling  stratification  on  the  northern  Skagit  Bay  tidal  flats.  The  authors  are 
Vera  Pavel,  B.  Raubenheimer,  and  Steve  Elgar.  The  appendix  is  part  of  this  manuscript. 


3.1  Introduction 

Stratification  has  been  studied  in  estuaries  and  coastal  river  plumes,  but  rarely  on  tidal 
flats  where  the  ratio  of  tidal  amplitude  to  depth  is  large  and  the  bed  is  dry  at  low  tide 
(effectively  “resetting”  the  system).  In  addition,  the  Skagit  tidal  flats  are  short  and  wide, 
have  freshwater  discharge  at  both  the  north  and  south  ends,  and  have  tidal  propagation 
parallel  to  the  bathymetry  contours,  unlike  many  estuaries  that  are  relatively  long  and 
narrow.  In  this  chapter,  field  observations  of  flows,  water  density,  and  water  levels  on  the 
Skagit  Bay  tidal  flats  will  be  used  to  evaluate  the  processes  creating  and  destroying 
stratification  during  tides  with  small  and  large  diurnal  inequalities,  The  relative  strength 
of  tidal  and  density-driven  processes  will  be  examined.  The  stratification  and  the 
processes  that  control  it  affect  turbulence  and  particle  transport,  and  thus  the  results 
presented  here  will  improve  our  understanding  of  the  hydrodynamic  and  ecological 
development  of  tidal  flats. 

3.2  Measurements 

Measurements  of  water  level,  currents,  and  water  density  were  collected  between  1 8  and 
31  August  2009  at  5  locations  (symbols  in  Fig.  2-1)  perpendicular  to  and  along 
(northwest  to  southeast)  the  2.5-m  depth  contour  on  the  tidal  flats.  Water  density  was 
estimated  from  measurements  with  induction-type  conductivity-temperature-depth  (CTD) 
sensors.  Near-bed  density  was  measured  with  a  fixed  CT  located  0.4  m  above  the  bed. 
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Near-surface  density  was  measured  with  a  CT  and  CTD  mounted  on  a  pole  at  distances 
of  0.2  and  0.7  m  below  a  surface  float,  respectively.  The  tilt  of  the  pole,  and  thus  the 
depth  of  the  upper  CT  sensor,  was  estimated  from  the  along-pole  distances  and  the  depth 
measured  by  the  lower  CTD.  Laboratory  tests  showed  that  errors  in  temperature  and 
salinity  are  less  than  about  0.1  °C  and  0.1  PSU,  respectively.  The  resulting  density 
accuracy  is  about  ±0.1  kg/m  and  the  depth  accuracy  is  about  ±0.01  m.  Bottom  pressure 
was  measured  at  4  Hz  with  pressure  sensors  buried  about  0.1  m  below  the  bed  level. 
Atmospheric  pressure  was  measured  at  4  Hz  near  La  Conner,  WA.  Nearbed  flows  were 
measured  about  0.1  m  above  the  bed  with  acoustic  Doppler  velocimeters  (ADVs) 
(accuracy  about  ±0.01  m/s)  that  collected  3072  s  of  data  at  2  Hz  starting  at  the  beginning 
of  each  hour.  Flow  profiles  were  measured  at  about  2  Hz  in  0.25-m  bins  from  about  0.4  m 
above  the  bed  to  one  to  two  bin-sizes  from  the  water  surface  with  2  MHz  upward-facing 
current  profilers  (accuracy  about  ±0.03  m/s  for  1-min  averages).  Instrument  locations 
were  surveyed  with  post-processed  differential  GPS  (accuracy  about  0.03  m).  Instruments 
were  separated  in  the  cross-  and  alongshore  by  about  600  and  1600  m,  respectively. 
Tidally  averaged  salinity  fields  from  model  simulations  suggest  that  alongshore  spatial 
scales  of  salinity  gradients  are  ~2-10  times  as  long  as  cross-shore  scales. 

Density  measurements  were  averaged  over  512-s  periods.  At  the  central  location,  the 
upper  floating  CT  failed,  and  the  data  were  replaced  by  an  average  of  the  4  other  upper 
CTs.  Root-mean-square  differences  between  the  512-s  averaged  densities  from  individual 
CT  measurements  and  the  average  of  the  values  from  all  sensors  at  the  same  elevation 
were  about  2  kg/m  ,  or  about  25%  of  the  surface-to-bottom  density  difference. 

Differences  between  the  sensor-averaged  density  and  the  measurement  at  the  central 
location  for  the  lower  floating  CTD  sensor  also  were  about  2  kg/m  .  Results  were  similar 
for  averages  including  data  from  only  the  alongshore  (on  the  same  depth  contour)  or  only 
the  cross-shore  sensors.  Atmospheric  pressure  was  removed  from  measured  bottom 
pressures.  These  adjusted  pressures  were  averaged  over  512  s,  and  used  to  estimate  water 
levels  assuming  hydrostatic  pressure  and  using  water  density  measured  by  the  CT 
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sensors.  Pressure  drifts  were  then  removed  by  subtracting  a  quadratic  fit  to  low-tide  data, 
when  the  flat  was  dry  and  the  water  depth  should  have  been  negligible.  AD  Vs  were 
assumed  to  be  fouled  or  out  of  the  water  and  data  were  discarded  when  the  signal  strength 
was  low  or  flows  were  noisy  (defined  as  times  when  the  root-mean-square  velocity 
fluctuations  were  more  than  twice  the  fluctuations  expected  from  applying  linear  theory 
to  the  pressure  signal).  Flows  measured  with  ADVs  were  averaged  over  512  s,  ignoring 
any  points  where  the  data  were  discarded.  The  estimated  water  levels  were  used  to 
determine  which  profiler  bins  were  actually  giving  measurements  of  the  water  column. 
Current  profiles  based  on  profiler  data  were  averaged  over  600-s  periods,  then 
interpolated  in  time  and  output  every  512  s  to  combine  with  ADV  data. 


Time  (day  in  August) 

Figure  3-1 :  (a  and  b)  Water  density  (color  contours)  as  a  function  of  depth  (thick  black  curve) 
and  time,  and  (c  and  d)  <I>  and  (e  and  f)  3<b/dt  vs.  time  for  selected  times  of  type  1  (a,  c,  and  e) 
and  type  2  (b,  d,  and  f)  tides. 

Density  is  dominated  by  salinity,  which  ranges  from  fresh  river  water  (density  ~1000 
kg/m  )  to  Puget  Sound  water  with  salinity  28  PSU  (density  ~1  020  kg/m  ).  Cross-shore 
and  alongshore  density  gradients,  based  on  differences  between  sensors,  had  median 
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values  of  about  5xl0'3  and  5xl0'4  kg/m4,  respectively.  Maximum  water  depths  ranged 
from  about  4  m  at  the  most  offshore  sensor  to  about  2  m  at  the  most  onshore  sensor. 
Cross-shore  and  alongshore  velocity  ranged  between  about  ±0.5  m/s. 


Time  (hours) 

Figure  3-2:  Phase-averaged  (a  and  b)  cross-  and  (c  and  d)  alongshore  velocity  (color  contours)  as 
a  function  of  depth  (thick  black  curve)  and  time  for  type  1  (a  and  c)  and  type  2  (b  and  d)  tides 


Similar  to  prior  observations  of  salt-wedge  estuaries  (Ralston  et  al.,  2010;  Giddings  et  al., 
2011)  and  numerical  simulations  of  Skagit  Bay  (Yang  and  Khangaonkar,  2009),  the  water 
column  is  fresh  at  the  beginning  of  flood  tide,  just  after  the  tidal  flat  is  submerged  (Fig.  3- 
1).  The  salinity  front  moves  onshore  during  flood,  with  the  water  column  becoming 
increasingly  saline.  During  ebb,  the  salinity  front  moves  offshore  and  the  water  freshens 
while  draining  off  the  tidal  flats. 
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During  type  1  tides,  cross-  and  alongshore  velocities  (Figs.  3-2a  and  c)  have  similar 
magnitudes  (see  also  Webster  et  al.,  2012).  Cross-shore  flows  are  relatively  weak  during 
type  2  tides  (Fig.  3-2b),  which  correspond  to  neap  tides  during  this  time  period,  but 
alongshore  flows  are  large  (Fig.  3 -2d),  resulting  in  similar  total  flow  magnitudes  during 
both  types  of  tides. 


3.3  Theory  and  Processing 


3.3.1  Theory 

The  stratification  is  quantified  using  the  potential  energy  anomaly  (O),  the  amount  of 
energy  per  unit  volume  required  to  homogenize  the  water  column: 

®  =  -j-fp'zdz,  (1) 

where  g  is  the  gravitational  acceleration,  h  and  r|  are  the  mean  depth  and  the  elevation  of 
the  surface  above  the  mean,  D=h+r\  is  the  total  water  depth,  p  =  p  +  p'  is  the  density, 
where  p  and  p'  are  the  depth-mean  and  residual  values,  and  z  is  the  vertical  coordinate, 
which  is  zero  at  the  mean  surface  and  positive  upward. 


The  temporal  evolution  of  <t>  (Fig.  3-1)  is  given  by  (Burchard  and  Hofmeister  2008): 
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where  t  is  time,  Vh  indicates  the  horizontal  components  of  the  gradient  operator,  u  is  the 
vector  horizontal  velocity  where  the  overbar  denotes  the  depth-average  and  the  prime  is 
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the  deviation  from  the  depth-average,  po  is  a  reference  density,  P§  and  are  surface  and 
bottom  buoyancy  flux,  respectively,  Q  represents  inner  sinks  and  sources,  Kh  is  the 
horizontal  eddy  diffusivity  and  w  is  the  deviation  of  the  vertical  velocity  (w)  from  a 
linear  profile: 


vi>  =  w-  —  +  u- 
\  dt 


z  +  h  _  )]  -  z 
- u-VJi- - 


The  eddy  diffusivity  Kv  is  estimated  as  (Munk  and  Anderson,  1948;  Nepf  and  Geyer, 
1996;  Burchard  and  Hofmeister,  2008;  Becker  et  al.,  2009): 

tfv-tf0(l  +  3.33^)~*  (4) 

where  R,  is  the  bulk  Richardson  number  (Byun  and  Wang,  2005;  Stacey  and  Ralston, 
2005)  and  K0  is  the  estimated  eddy  diffusivity  for  an  unstratified  water  column  (see 
Appendix  for  further  details  on  the  mixing  parameterization).  The  last  three  terms  on  the 
right  side  of  equation  (2)  are  expected  to  be  small  compared  with  the  other  terms  and  are 
neglected.  To  examine  the  effects  of  changing  depth,  the  first  term  is  separated  into  an 
advection  and  depth  change  term: 

-Vh(u<S>)  =  -uVh<S>-<S>Vhu  (5) 

By  continuity,  the  horizontal  gradient  of  the  mean  velocity  can  be  expressed  in  terms  of 
depth  changes  as: 

V  ,u  =  - — f  —  +  0  •  V  ,d\  (6) 


D\dt 


Finally,  rearranging  the  integrals  as: 


- z 

2 


Xdz  =  f  zx'dz 


where  x  is  an  arbitrary  depth-dependent  function,  the  overbar  denotes  the  depth-average, 
and  the  prime  denotes  the  deviation  from  the  depth-average,  results  in: 
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where  term  A  is  direct  advection  of  stratification,  DC  is  the  effect  of  depth  changes  in 
both  time  and  space,  DS  is  depth-mean  straining,  NS  is  non-mean  straining  caused  by 
shear  acting  on  a  non-depth-uniform  density  gradient,  M  is  mixing,  and  VA  is  vertical 
advection  shifting  the  isopycnals  up  and  down, 


3.3.2  Processing 

To  evaluate  the  integrals  in  equation  (8),  density  and  velocity  profiles  were  linearly 
interpolated  onto  a  0.1 -m  vertical  grid.  Profiles  were  extended  to  the  surface  and  bed 
assuming  constant  values  given  by  the  highest  and  lowest  measurement,  respectively. 
Results  are  similar  to  assuming  zero  flux  at  the  bed  and  surface  and  using  a  linear 
extrapolation  of  the  vertical  gradient.  Errors  in  O  owing  to  the  vertically  sparse  density 
measurements  are  of  order  10%,  based  on  differences  between  this  method  and  fitting  to 
a  two-layer  model.  Near-bottom  density  was  interpolated  along  the  sloping  bed.  The 
vertical  structure  of  cross-shore  density  gradients  at  elevations  between  the  bed  levels  at 
any  two  locations  was  calculated  using  the  density  profile  at  the  deeper  location  and  the 
corresponding  densities  (same  elevations)  estimated  along  the  sloping  bed  (Fortunato  and 
Baptista,  1996).  Horizontal  gradients  were  evaluated  using  upstream  (determined  from 
the  depth-averaged  velocity  at  the  central  location)  differences. 

Profiler  measurements  of  horizontal  flows  coupled  with  a  mass  balance  suggest  that 
vertical  velocities  (w)  were  smaller  than  the  resolution  of  the  instrument  (0.001  m/s),  so 
the  vertical  advection  term  could  not  be  calculated  from  data.  Previous  results  suggest 
that  this  term  is  unlikely  to  be  large  except  in  a  small  region  near  the  density  front  (Nepf 
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and  Geyer,  1996;  Burchard  and  Hofmeister,  2008;  de  Boer  et  al.,  2008;  Hofmeister  et  al., 
2009;  Marques  et  al.,  2010),  and  thus  it  is  neglected  here.  However,  vertical  advection 
may  contribute  to  the  errors  in  the  estimated  stratification  balance. 


The  water  depth,  d>,  and  all  terms  in  equation  (8)  were  smoothed  using  a  7200-s  running 
average.  The  averaging  period  was  chosen  to  be  shorter  than  %-tidal  period  to  resolve 
tidal  fluctuations,  but  longer  than  the  advective  timescale  between  sensor  locations  (about 
3600  s  and  4800  s  in  the  cross-  and  alongshore,  respectively).  This  temporal  averaging 
and  smoothing  reduces  errors  owing  to  unresolved  small-scale  spatial  variability.  O  was 
smoothed  before  difference  calculations  were  performed  and  the  results  were  smoothed 
again  after  computing  dO/dt  and  the  advection  term  to  remove  jitter  artifacts  that  result 
from  numerical  differencing  of  noisy  data.  All  terms  were  phase  averaged  over  24-h  long 
(diurnal)  cycles  for  type  1  (5  cycles)  and  type  2  (4  cycles)  tides.  The  distinction  between 
the  tide  types  1  and  2  was  made  by  comparing  the  difference  in  water  level  of  the  two 
low  tides  to  the  total  tidal  range  for  that  diurnal  cycle.  If  the  difference  was  less  than  1/4 
the  tidal  range  it  was  considered  a  type  1  tide.  If  the  difference  was  greater  than  1/3  the 
tidal  range  it  was  considered  a  type  2  tide. 


3.4  Results 

3.4.1  The  Potential  Energy  Anomaly  Balance 

The  ensemble-averaged  diurnal  variations  of  dO/dt  (Fig.  3-3),  including  the  timing  of  the 
maxima  and  minima,  are  consistent  with  the  stratification  balance  (equation  (8)).  The 
good  agreement  between  the  left  hand  side  of  the  equation  (Fig.  3-3c  and  d,  black  curve), 
which  was  calculated  directly  from  the  central  site’s  O  (Fig.  3-3a  and  b),  and  the  right 
side  of  the  balance  (Fig.  3-3c  and  d,  purple  dashed  curve)  suggests  that  the  sum  of  the 
neglected  terms  (vertical  advection,  additional  sources  of  mixing,  and  surface  or  bed 
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buoyancy  fluxes)  is  not  large  compared  with  the  retained  terms,  possibly  because  the 
neglected  terms  cancel,  as  in  numerical  simulations  of  the  Patos  Lagoon  ROFI  (Marques 
et  al.,  2010). 
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Figure  3-3:  Phase-averaged  <I>  (a  and  b)  and  dd>/dt  (c  and  d)  based  on  the  observed  density 
profiles  (solid  black  curves)  and  on  estimations  of  stratification-related  processes  (right-hand  side 
of  equation  (1),  dotted  purple  curves)  versus  time  for  (a  and  c)  type  1  and  (b  and  d)  type  2  tides. 
Squared  correlations  between  the  direct  and  process  estimates  are  about  0.7.  The  shaded  area 
shows  the  relative  water  depth  over  the  tidal  cycle. 


During  the  type  1  and  2  strong  flood  tides  the  thin  layer  of  water  initially  submerging  the 
flats  is  slightly  stratified  (Figs.  3-3a  and  b),  but  becomes  increasingly  well-mixed  (Figs. 
3-3c  and  d,  dO/dt  <  0.  In  contrast  to  many  partially  mixed  estuaries  and  ROFIs  in  which 
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depth-mean  straining  dominates  and  maximum  stratification  occurs  during  late  ebb  or 
low  water  (Nepf  and  Geyer,  1996;  Rippeth  et  al.,  2001;  Burchard  and  Hofmeister,  2008), 
but  similar  to  observations  in  strongly-forced  salt- wedge  estuaries  (Ralston  et  al.,  2010a; 
Giddings  et  al.,  2011),  maximum  stratification  occurs  at  about  mid  ebb  tide  (Fig.  3-3c 
time  about  9  and  21  h,  where  dO/9t  changes  from  positive  to  negative).  During  type  2 
tides,  stratification  increases  during  the  beginning  of  the  weak  flood  (Fig.  3-3d,  time  from 
about  13  to  16  h,  dO/dt  >  0,  and  Fig.  3-4b).  Although,  temporal  changes  in  stratification 
vary  during  individual  type  2  tidal  cycles,  in  3  of  4  cases  O  is  larger  at  the  end  of  the 
weak  flood  than  at  the  beginning  of  the  weak  ebb  (Fig.  3-4b). 

In  prior  studies  of  stratification  for  mixed  tides  (Ralston  and  Stacey,  2005b;  Wang  et  al., 
2009;  Giddings  et  al.,  201 1),  the  stratification  decreases  during  the  weak  flood,  briefly 
increases  during  the  beginning  of  the  following  strong  ebb,  then  decreases  again,  which 
appears  to  be  a  result  of  straining  on  the  weak  flood  and  beginning  of  ebb,  that  is  then 
combined  with  mixing,  and  possibly  advection,  later  on  the  ebb.  The  measurements  in  the 
Skagit  for  similar  time  periods  (Fig.  3 -4b  time  ~16  to  22  h)  show  high  variability,  but  on 
average,  the  stratification  decreases  from  about  mid  weak  flood  through  the  second  high 
tide  and  the  following  strong  ebb.  Although  the  same  processes  are  likely  acting,  this 
difference  suggests  that  the  advection  process  may  be  more  influential  in  the  Skagit  for 
this  tidal  phase. 

The  individual  terms  in  the  balance  (Fig.  3-5)  indicate  which  processes  dominate  the 
stratification.  The  fresh  water  draining  off  the  flats  during  the  strong  ebbs  remains  partly 
stratified  (Figs.  3-lc  and  d).  Thus,  on  the  strong  floods,  advection  and  depth-mean 
straining  (solid  blue  and  dashed  purple  curves  in  Fig.  3 -5a  for  time  about  2  to  6  and  12  to 
16  h,  and  Fig.  3 -5b  for  time  about  2  to  6  h)  often  are  negative  as  the  thin,  mostly  fresh 
tongue  that  initially  covers  the  flats  is  replaced  by  water  from  offshore  that  is  increasingly 
well-mixed  and  saline  (Figs.  2-3  b  and  d).  Advection  usually  opposes  the  positive  depth- 
mean  straining  on  the  ebb  (Fig.  3-5  time  about  8  to  12  and  18  to  22  h),  indicating  that  the 
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fresh  water  that  is  trapped  near  the  shore  at  high  tide  (Fig.  2-3  d)  is  less  stratified  than  the 
offshore  water,  similar  to  many  estuaries  (Burchard  and  Hofmeister,  2008;  Ralston  et  ah, 
2010a;  Giddings  et  ah,  2011)  where  river  discharge  occurs  onshore,  but  in  contrast  to 
models  of  ROFIs  (de  Boer  et  al.,  2008)  where  the  freshwater  source  is  not  directly 
onshore.  The  negative  advection  also  suggests  that  stratified  water  is  removed  from  the 
upper  flats  and  transported  to  the  offshore  flats  during  the  strong  ebb,  where  it  would  be 
transported  northward  towards  Deception  Pass  (Fig  2-3e).  Comparison  of  the  terms 
during  type  1  tides  with  those  during  type  2  tides  suggests  that  large  negative  values  of 
dO/dt  occur  on  strong  (but  not  weak)  ebbs  at  least  partly  because  advection  of  well-mixed 
water  is  larger  on  strong  ebbs  (compare  the  maximum  negative  values  of  the  dark  blue 
curve  in  Fig.  3-5a  with  those  in  Fig.  3-5b,  e.g.,  near  time  is  10  h).  The  larger  advection  on 
strong  ebbs  than  weak  ebbs  is  a  primarily  a  result  of  the  difference  in  the  strength  of  the 
tidal  flows. 
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Figure  3-4:  Comparison  of  individual  tidal  cycles  (thin  curves)  to  the  average  (thick  black  curve) 
for  O  (a  and  b)  and  dQ/dt  (c  and  d)  for  type  1  (a  and  c)  and  type  2  (b  and  d)  tides. 
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Near  the  end  of  the  type  2  weak  ebb,  advection  is  small  and  switches  from  negative  to 
positive  (dark  blue  curve  in  Fig.  3-5b  near  time  about  12  h),  at  least  partly  owing  to 
positive  alongshore  advection  (Fig.  3 -6b)  caused  by  strong  southeasterly  alongshore 
flows  (Figs.  2-3e  and  3-2d)  transporting  stratified  water  alongshore  from  the  north  fork  of 
the  Skagit  River  towards  the  instrument  location  and  the  central  flats  (Fig  3-7b).  The 
water  on  the  offshore  flats  remains  stratified  following  the  weak  ebb  (Fig.  3-7d),  and  thus 
cross-shore  advection  is  positive  throughout  the  weak  flood.  Cross-shore  depth-mean 
straining  also  remains  positive  during  the  small  low  and  until  about  mid  weak  flood  (Fig. 
3-5b)  owing  to  vertically  sheared  currents  with  offshore  flow  near  the  surface  (Fig.  3-2b). 

Mixing  (Fig.  3-5,  green  dashed  curve),  which  always  is  negative  (reducing  stratification), 
is  largest  when  internal  or  bottom  shear  is  large  and  when  the  water  is  not  well  mixed 
already.  During  type  1  tides,  mixing  is  strongest  during  the  ebbs,  which  is  consistent  with 
the  combined  interfacial-  and  bottom-generated  mixing  observed  in  strongly  forced  salt- 
wedge  estuaries  (Ralston  et  al.,  2010b;  Giddings  et  al.,  2011;  Wang  et  al.,  2011).  Mixing 
is  weak  during  floods  despite  strong  near-bed  flows  primarily  because  there  is  a  lack  of 
stratification  to  be  mixed.  However,  complete  mixing  does  not  seem  to  occur  and  the 
minimum  (geometrically)  phase-averaged  bulk  Richardson  number  is  about  0.2, 
occurring  after  the  strong  flood  of  type  1  tides.  During  type  2  tides,  mixing  is  relatively 
large  during  the  small  low  tide  (dashed  green  curve  in  Fig.  3-5b  for  time  about  12  h), 
similar  to  observations  and  model  predictions  in  the  Snohomish  estuary  (Wang  et  al., 
2009;  Giddings  et  al.,  2011).  Near-bed  flows  are  weak,  while  mid-water  column  shear  is 
strong  (Figs.  3-2b  and  d),  suggesting  the  mixing  on  type  2  tides  primarily  is  caused  by 
internal  shear,  especially  for  the  more  highly  vertically  sheared  alongshore  currents  (Fig. 
3-2d).  The  bulk  Richardson  number  at  the  higher  low  water  is  about  0.5,  which  is  the 
lowest  value  during  phase-averaged  type  2  tides,  equal  to  that  occurring  at  the  end  of  the 
type  2  strong  flood.  Similar  to  other  systems  with  strong  stratification  (Burchard  and 
Hofmeister,  2008;  de  Boer  et  al.,  2008),  mixing  usually  is  smaller  than  advection  and 
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depth-mean  straining  (Fig.  3-5),  but  here  mixing  remains  significant  owing  to  the  shallow 
depths. 

The  non-mean  straining  is  small  during  type  1  tides  (Fig.  3-5a,  solid  light  blue  curve). 
However,  in  contrast  to  numerical  simulations  showing  that  non-mean  straining  typically 
is  large  only  near  the  river  or  estuarine  mouth  (de  Boer  et  al.,  2008;  Marques  et  al.,  2010), 
here  it  is  similar  in  magnitude  to  the  other  terms  during  type  2  tides  (Fig.  3-5b),  and  it 
reduces  the  increase  in  stratification  during  the  small  low  tide  (Fig.  3-3d,  dO/dt  >  0)  and 
the  decrease  in  stratification  during  the  latter  half  of  the  weak  flood  (Fig.  3-3d,  dO/dt  < 

0).  The  prior  simulations  suggest  vertical  advection  may  balance  non-mean  straining,  and 
thus  the  discrepancies  between  the  direct-  and  process-based  estimates  of  dO/dt  (Fig.  3- 
3d)  during  the  small  low  tide  and  weak  flood  may  result  partly  from  neglecting  vertical 
advection.  Discrepancies  also  may  occur  because  the  depth  dependence  of  density  (which 
affects  non-mean  straining)  is  not  well-resolved  by  the  field  measurements. 

Similar  to  long,  narrow  estuaries  (Simpson  et  al.,  1990;  Nepf  and  Geyer,  1996;  Burchard 
and  Hofmeister,  2008;  Giddings  et  al.,  2011),  the  cross-shore  advection  (Fig.  3-6a)  and 
depth-mean  straining  (Fig.  3-6c)  during  type  1  tides  are  5  and  15  times  larger, 
respectively,  than  the  alongshore  components.  However,  the  cross-  and  alongshore 
velocities  (Figs.  3 -2a  and  c)  have  similar  magnitudes,  consistent  with  the  dominance  of 
cross-shore  processes  being  a  result  of  larger  cross-shore  density  gradients  (see  also  Figs. 
2-3  and  2-4).  The  depth-mean  (Fig.  3-6d)  and  non-mean  (Fig.  3-6f)  straining  also  are 
dominated  by  the  cross-shore  component  during  type  2  tides,  but  the  alongshore 
component  of  the  advection  during  type  2  tides  has  about  half  the  magnitude  of  its  cross¬ 
shore  component  (Fig.  3-6b).  The  largest  magnitudes  of  alongshore  advection  for  type  2  t 
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Figure  3-5  Phase-averaged  values  of  terms  on  the  right  side  of  equation  (2)  versus  time  for  (a) 
type  1  and  (b)  type  2  tides.  The  shaded  area  shows  the  relative  water  depth  over  a  tidal  cycle. 
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ides  result  from  large  alongshore  stratification  variations  (Fig.  3 -7b  time  from  10  to  16  h) 
and  alongshore-dominated  flows  (compare  Figs.  3-2d  with  3-2b  for  time  from  10  to  16  h) 
during  the  end  of  the  weak  ebb  and  the  beginning  of  the  weak  flood  owing  to  the  large- 
scale  tidal  circulation  (Yang  and  Khangaonkar,  2009).  Owing  to  the  difference  in  tidal 
range  at  Deception  Pass  and  Saratoga  Passage,  there  is  a  flow  to  the  northwest  at  high 
tides  and  a  flow  to  the  southeast  at  low  tides.  The  high  tide  effect  is  not  seen  in  the 
stratification  balance  of  type  1  tides,  and  is  smaller  than  the  low  tide  effect  during  type  2 
tides  because  alongshore  density  and  stratification  gradients  are  small  during  the  high 
tide,  a  result  of  the  measurement  location  being  well  inside  the  unstratified  salt  wedge. 
The  low  tide  effect  is  not  seen  on  type  1  tides  because  there  is  no  water  on  the  flats  at  the 
measurement  location. 


3.4.2  Spatial  Variability 

As  in  prior  studies  (Burchard  and  Hofmeister,  2008;  de  Boer  et  al.,  2008;  Marques  et  ah, 
2010),  the  potential  energy  anomaly  and  the  dominant  processes  controlling  the 
stratification  can  be  spatially  variable  (Fig.  3-7).  Here,  the  strength  of  the  stratification 
typically  increases  towards  the  northwest  (towards  the  mouth  of  the  north  fork,  Figs.  3-7a 
and  b)  and  offshore  (Figs.  3-7c  and  d).  For  both  type  1  and  2  tides,  the  local  maxima  in 
stratification  occur  earlier  at  onshore  locations  (red  dotted  curves  in  Figs.  3-7c  and  d) 
than  at  offshore  locations  (purple  solid  curves),  possibly  owing  to  the  earlier  passage  of 
the  density  front. 

During  type  1  tides,  temporal  changes  in  are  similar  at  all  sensors  (Figs.  3-7a  and  c),  in 
contrast  to  larger,  more  spatially  variable  estuaries  and  ROFIs  (de  Boer  et  ah,  2008; 
Burchard  and  Hofmeister,  2008;  Wang  et  ah,  2009;  Marques  et  ah,  2010,  201 1;  Ralston 
et  ah,  2010a;  Giddings  et  ah,  2011).  The  front-induced  stratification  maximum  on  the 
flood  following  the  smaller  low  becomes  stronger  at  the  northwest  (closer  to  the  north 
fork,  Fig.  3-7a  for  time  about  20  to  24  h)  and  offshore  (Fig.  3-7c  for  time  about  20  to  24 
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h)  locations,  possibly  because  the  salt  wedge  becomes  stronger  closer  to  the  river  mouth 
and  offshore. 


Time  (hours) 


Figure  3-6:  Phase-averaged  cross-  (solid  curves)  and  alongshore  (dotted  curves)  components  of 
(a  and  b)  advection,  (c  and  d)  depth-mean  straining,  and  (e  and  f)  non-mean  straining  for  (a,  c, 
and  e)  type  1  and  (b,  d,  and  f)  type  2  tides  versus  time.  The  shaded  area  (g  and  h)  shows  the 
relative  water  depth  over  a  tidal  cycle. 


The  maximum  stratification  is  larger  during  type  2  tides  than  during  type  1  tides  at  most 
locations  (Fig.  3-7).  Consistent  with  prior  observations  in  shallow  channels  and  in 
strongly  forced  salt-wedge  estuaries  (Ralston  and  Stacey,  2005b;  Wang  et  al.,  2009; 
Ralston  et  al.,  2010a;  Giddings  et  al.,  2011),  the  water  column  remains  stratified  during 
the  type  2  small  low  tide.  The  density-driven  circulation  dominates  during  these  weak 
tidal  flows,  except  at  the  most  onshore  location,  which  may  have  uniform  density  at  that 
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time  (Fig.  3-7d).  The  stratification  increases  during  the  weak  flood  at  most  locations 
(Figs.  3-7b  and  d),  presumably  owing  to  the  same  advection  of  stratified  water  (Fig.  3-6b) 
and  depth-mean  straining  (Fig.  3-6d)  observed  at  the  central  location.  Except  at  the 
southeastern  and  onshore  locations,  the  water  column  remains  partly  stratified  until  water 
drains  off  the  flat. 


Time  (hours) 


Figure  3-7:  Phase-averaged  potential  energy  anomaly  for  (a  and  b)  alongshore  and  (c  and  d) 
cross-shore  sensor  lines  for  (a  and  c)  type  1  and  (b  and  d)  type  2  tides,  and  (e  and  f)  water  depth 
versus  time. 


The  spatially  uniform  temporal  evolution  of  stratification  during  type  1  tides  suggests  that 
the  processes  affecting  stratification  are  similar  across  the  flats.  However,  the  spatial 
variability  during  type  2  tides  suggests  that  different  processes  may  dominate  depending 
on  location,  similar  to  numerical  model  simulations  over  large  regions  of  salt- wedge 
estuaries  (Wang  et  al.,  2009;  Ralston  et  al.,  2010a)  and  ROFIs  (Marques  et  al.,  2010).  The 
difference  between  type  1  and  type  2  stratification  is  at  least  partly  owing  to  the  relative 
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importance  of  tidal  barotropic  circulation,  which  affects  the  advection  and  mixing.  In 
particular,  barotropic  forcing  is  larger  than  density-driven  estuarine-style  circulation 
during  the  strong  floods  and  ebbs.  In  contrast,  during  the  type  2  weak  flood  and  ebb,  the 
estuarine  circulation  increases  stratification  with  little  modulation  from  the  weakened 
barotropic  forcing.  During  this  period,  spatial  variations  in  stratification  are  relatively 
large  because,  unlike  the  externally  forced  tides,  the  strength  of  estuarine  circulation  is 
sensitive  to  local  density  gradients  and  water  depth. 


3.4.3  Robustness  of  results 


Type  1  tides  Type  2  tides 

Min 

Mean 

Max 

Min 

Mean 

Max 

0.75 

0.83 

0.86 

0.67 

0.84 

0.95 

d<D/d/ 

0.76 

0.84 

0.95 

0.49 

0.66 

0.83 

Advection 

0.53 

0.71 

0.80 

0.34 

0.68 

0.83 

Depth-mean  straining 

0.77 

0.85 

0.93 

0.66 

0.74 

0.86 

Non-mean  straining 

0.37 

0.58 

0.72 

0.34 

0.68 

0.93 

Mixing 

0.40 

0.59 

0.89 

0.62 

0.71 

0.76 

Depth  Change 

0.48 

0.71 

0.88 

0.57 

0.80 

0.96 

Table  3-1:  Correlation  coefficients  as  r2.  Minimum,  mean,  and  maximum  are  given  for  individua 
tidal  cycles  compared  against  the  phase  averaged  results  given  in  Figs.  3-3  and  3-5  . 


Correlations  of  d>,  directly  calculated  dO/dt,  and  each  of  the  terms  shown  in  Fig.  3-5 
were  performed  between  individual  tidal  cycles  (Fig.  3-4)  and  the  phase-averaged  results 
(Table  3-1).  Most  of  the  values  of  r  were  above  0.6,  and  many  were  above  0.8,  which 
shows  that  the  phase-averaging  method  preserves  most  of  the  variability  of  these 
processes.  Depth-mean  straining,  O,  and  dO/dt  had  particularly  high  correlations, 
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indicating  that  tidal  variations  of  these  quantities  are  more  stable  between  cycles.  For 
type  1  tides  correlations  are  similar  for  dO/dt  and  O,  while  for  type  2  tides  the  O 
correlation  is  higher.  For  type  1  tides,  both  O  and  dO/dt  have  2  or  3  strong  maxima  and 
minima  associated  with  the  strong  floods  and  ebbs  (Fig.  3-4a  and  c).  In  contrast,  during 
type  2  tides,  O  rises  on  the  weak  ebb  and  remains  high  until  the  strong  ebb  (Fig.  3-4b). 
Small  fluctuations  during  this  period  cause  oscillations  in  dO /dt  that  are  different  in 
magnitude  and  phase  between  the  individual  tidal  cycles  (Fig.  3-4d). 

3.5  Summary  and  Conclusions 

The  key  features  of  the  different  tidal  stages  are  listed  below. 

Processes  occurring  on  strong  floods  during  type  1  and  type  2  tides  (Fig.  3-8)  include: 

•  Water  moves  onshore. 

•  O  is  positive  at  the  start  of  the  flood  as  the  salt  wedge  encroaches  (initially  near 
the  seafloor)  on  the  riverine  water  remaining  on  the  flats. 

•  The  velocity  shear  moves  the  denser  offshore  water  onshore  faster  at  the  surface 
than  near  the  bed,  reducing  O  through  (negative)  depth-mean  straining. 

•  The  offshore,  previously  unstratified,  water  is  advected  onshore. 
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Strong  Flood 


Side  view 


Figure  3-8:  Cross-section  (a)  and  plan  view  (b)  diagrams  of  the  flats  during  the  strong  flood.  The 
white  arrows  (a)  indicate  the  cross-shore  flow  profile,  and  colored  arrows  (b)  indicate  the 
stratification  processes  with  onshore-  and  offshore-directed  arrows  indicating  increasing  and 
decreasing  stratification,  respectively.  The  length  of  each  arrow  indicates  the  relative  strength  of 
the  flow  or  process.  Blue  shading  represents  water,  with  lighter  colors  representing  fresher  water. 
The  cross-hatching  indicates  the  stratified  region. 


Processes  occurring  on  type  1  strong  ebbs  (Fig.  3-9)  are: 

•  Water  moves  offshore. 

•  O  initially  is  small,  then  increases  strongly,  and  then  decreases. 

•  During  the  increase  in  O,  the  velocity  shear  moves  the  fresher  water  from  onshore 
to  offshore  faster  at  the  surface  than  near  the  bed,  increasing  O  through  depth- 
mean  straining 

•  O  created  by  the  depth-mean  straining  tends  to  be  larger  in  deeper  water,  so  as  the 
ebb  progresses,  the  low-O  onshore  water  is  advected  offshore  by  the  mean  flow, 
resulting  in  negative  values  of  advection. 

•  Mixing  also  becomes  relevant  on  the  ebb  in  the  presence  of  strong  stratification 
and  relatively  strong  flows. 
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Strong  Ebb 


Side  view 


Figure  3-9:  Cross-section  (a)  and  plan  view  (b)  diagrams  of  the  flats  during  the  type  1  strong  ebb. 
The  white  arrows  (a)  indicate  the  cross-shore  flow  profile,  and  colored  arrows  (b)  indicate  the 
stratification  processes  with  onshore-  and  offshore-directed  arrows  indicating  increasing  and 
decreasing  stratification,  respectively.  The  length  of  each  arrow  indicates  the  relative  strength  of 
the  flow  or  process.  Blue  shading  represents  water  with  lighter  colors  representing  fresher  water. 
Stratification  is  present  throughout  the  flats  and  is  not  noted. 

Processes  occurring  during  the  extended  inundation  (weak  ebb,  higher  low  water,  and 
weak  flood)  of  type  2  tides  (Fig.  3-10): 

•  On  average,  O  increases  quickly  on  the  weak  ebb  and  continues  to  increase 
slowly  through  higher  low  water  and  the  beginning  of  weak  flood,  before 
decreasing  at  the  end  of  the  weak  flood. 

•  The  tidal  velocity  is  relatively  weak,  and  density-driven  estuarine-style  circulation 
is  present,  which  results  in  positive  depth-mean  straining  throughout  this  time. 

•  O  is  higher  offshore  than  onshore,  resulting  in  cross-shore  advection  that  is 
negative  on  ebb  and  positive  on  flood. 

•  Mixing  peaks  at  low  water,  corresponding  to  the  strongest  flows  during  the  high- 
O  period. 

•  Owing  to  alongshore  differences  in  O,  but  not  depth-averaged  density,  strong 
alongshore  flows  with  strong  shear  result  in  alongshore  advection,  but  not 
alongshore  straining. 
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The  processes  on  the  strong  ebb  of  type  2  tides  are  the  same  as  those  on  type  1  strong 
ebbs,  except  there  is  no  delay  between  the  onset  of  straining  and  the  onset  of  negative 
advection  and  mixing  because  the  water  column  starts  stratified  instead  of  unstratified. 

Overall,  the  stratification  observed  on  the  shallow  and  wide  tidal  flats  near  the  north  fork 
of  the  Skagit  River,  where  flows  are  similar  in  both  across-and  along-isopycnal 
directions,  has  similarities  to  both  deeper-water  nearshore  regions  of  fresh  water 
influence  (ROFIs)  where  strong  flows  are  perpendicular  to  the  density  gradient  (Rippeth 
et  al.,  2001;  de  Boer  et  al.,  2008;  Marques  et  al.,  2010,  2011),  and  to  narrow  tidal-flat 
channels  (Ralston  and  Stacey,  2005a, b)  and  strongly  forced  salt-wedge  estuaries  (Wang 
et  al.,  2009;  Ralston  et  al.,  2010a, b;  Giddings  et  al.,  2011)  where  tidal  flows  and  density 
gradients  are  aligned. 


Extended  inundation 


Side  view 


Figure  3-10:  Cross-section  (a)  and  plan  view  (b)  diagrams  of  the  flats  during  the  type  2  extended 
inundation.  The  white  arrows  (a)  indicate  the  cross-shore  flow  profile,  and  colored  arrows  (b) 
indicate  the  stratification  processes  with  onshore-  and  offshore-directed  arrows  indicating 
increasing  and  decreasing  stratification,  respectively.  The  length  of  each  arrow  indicates  the 
relative  strength  of  the  flow  or  process.  Blue  shading  represents  water  with  lighter  colors 
representing  fresher  water.  The  cross-hatching  (b)  indicates  the  stratified  region,  with  larger 
squares  indicating  higher  stratification.. 
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The  water  on  the  tidal  flats  did  not  become  completely  mixed  during  the  strong  ebb,  and 
thus  the  initial,  partially  stratified  tongue  of  water  crossing  the  flats  on  strong  floods  may 
have  some  leftover  stratification  from  the  previous  tidal  cycle  as  well  as  from  the 
intrusion  of  incoming  Puget  Sound  water  into  the  residual  river  water.  As  in  prior  studies, 
the  water  becomes  increasingly  saline  during  the  strong  flood,  and  for  nearly  semi¬ 
diurnal  tides,  the  maximum  stratification  occurs  during  mid  ebb  tide.  Throughout  the  tidal 
cycle,  changes  in  stratification  result  primarily  from  cross-shore  tidal  straining  and 
advection. 

Stratification  is  stronger  during  nearly  diurnal  (type  2)  tides  than  during  nearly  semi¬ 
diurnal  (type  1)  tides.  In  contrast  to  prior  estuarine  studies,  but  similar  to  many  ROFIs, 
alongshore  flows  are  stronger  than  cross-shore  flows,  and  the  increased  stratification  on 
type  2  tides  is  owing  to  both  alongshore  advection  of  stratified  water  from  the  north  fork 
of  the  Skagit  River,  and  straining  effects  of  density-driven  cross-shore  shear  flows  that 
remain  offshore-directed  at  the  surface  until  mid  weak  flood,  when  stratification  is 
maximum  (Fig.  3-10).  The  non-mean  straining  during  the  extended  high  water  of  type  2 
tides  indicates  that  the  isopycnals  are  not  parallel  during  that  tidal  stage,  while  the 
negligible  non-mean  straining  during  the  type  1  tides  suggests  that  they  are  during  that 
type  of  tide. 

The  large  advection  terms  suggest  that  strong  fronts  are  moved  across  and  along  the  flats. 
Positive  straining  (indicating  generation  of  stratification)  usually  is  larger  than  negative 
mixing  (indicating  destruction  of  stratification),  suggesting  water  that  becomes  stratified 
near  the  measurement  location  is  either  mixed  elsewhere  on  the  flats,  or  exported  to  the 
surrounding  basins. 

Temporal  changes  in  stratification  are  similar  across  and  along  the  flats  in  this  region. 
However,  stratification  increases  offshore  and  alongshore  towards  the  north  fork  (the 
river  mouth  closest  to  the  instruments). 
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Overall,  type  1  tidal  processes  can  be  roughly  described  as  a  salt  wedge  with  a  nearly 
vertical  front  that  is  transported  onshore  on  the  flood  tide,  which  is  then  strained  so  it  is 
more  stratified  and  covers  a  larger  area  while  it  is  advected  offshore  on  the  ebb.  During 
type  2  tides,  the  salt  wedge  again  passes  onto  the  flats  during  the  strong  flood,  but  during 
the  prolonged  high  water  of  the  weak  ebb  and  flood,  additional  cross-shore  straining  and 
alongshore  advection  occur  with  minimal  mixing  (owing  to  weak  mean  flows),  resulting 
in  increased  stratification.  Runaway  stratification  in  the  measurement  location  is 
prevented  by  water  draining  off  the  flats  during  the  strong  ebb,  and  probably  offshore 
mixing. 
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4  Numerical  model  simulations  of  stratification  processes 

4.1  Introduction 

The  observations  (Chapter  3)  provide  information  about  the  processes  affecting 
stratification.  However,  the  results  are  applicable  only  to  the  northern  flats  and  to  a  period 
of  low  discharge.  Here,  the  spatial  dependence  of  the  processes  controlling  the  creation 
and  destruction  of  stratification  are  examined  using  the  Finite  Volume  Coastal  Ocean 
Model  (FVCOM). 

4.1.1  FVCOM  model 

Water  level  fluctuations,  flows,  water  density,  and  mixing  were  simulated  from  July  1  to 
15,  2009,  using  FVCOM,  a  three-dimensional  finite  volume,  unstructured  grid,  primitive 
equation  model  (Chen  et  al.,  2003).  FVCOM  has  been  used  with  success  in  modeling  the 
Baeksu  tidal  flats  (Kim  and  Cho,  2011)  and  compares  well  with  field  observations  in  the 
Merrimack  estuary  (Ralston  et  al.,  2010a)  as  well  as  having  been  used  to  investigate 
circulation  in  the  Skagit  (Ralston  et  al.,  2012;  Yang  and  Khangaonkar,  2009). 

4.1.2  Model  setup 

The  model  domain  encompasses  Skagit  Bay  and  a  portion  of  the  surrounding  waters  of 
Saratoga  passage  and  the  Strait  of  Juan  de  Fuca.  Cell  sizes  range  from  15  m  within  Skagit 
Bay  to  500  m  near  the  open  boundaries.  The  model  is  similar  to  that  used  by  Yang  and 
Khangaonkar  (2009),  but  with  an  expanded  domain,  higher  resolution  on  the  flats,  and 
updated  bathymetry.  Uniform  sigma  coordinates  were  used  in  the  vertical  and  “wet”  and 
“dry”  mesh  elements  were  used  to  account  for  the  moving  water  level.  Turbulence  was 
modeled  with  a  k-s  turbulence  closure  scheme.  The  density  in  Skagit  Bay  is  dominated 
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by  salinity  fluctuations,  and  thus  the  temperature  was  assumed  uniform  (18  °C). 
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Figure  4-1:  Model  bathymetry  and  instrument  locations  on  the  north  flats  (yellow,  July,  2009) 
and  south-central  flats  (red,  June,  2009).  The  black  line  segment  is  the  transect  examined  in 
section  4.2.2. 


The  model  was  forced  at  open  boundaries  with  predicted  tidal  fluctuations  and  observed 
subtidal  water  level  fluctuations  from  the  tidal  reference  stations  of  Seattle  (#9447130) 
and  Port  Townsend  (#9444900),  and  measured  river  discharge  at  the  USGS  station  at  Mt. 
Vernon,  WA  (#12200500).  Winds  were  forced  using  the  regional  Weather  Research  and 
Forecasting  model.  A  more  detailed  description  of  the  FVCOM  setup  for  the  Skagit  Bay 
tidal  flats  is  given  by  Ralston  et  al.  (2012). 
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4.1.3  Model  verification 


North  flats,  July  2009 

Min 

Mean 

Max 

r2 

SS 

r2 

SS 

r2 

SS 

Depth 

0.97 

0.40* 

0.99 

0.94 

1.00 

1.00 

Cross-shore  Velocity 

0.37 

0.22 

0.61 

0.53 

0.78 

0.68 

Alongshore  Velocity 

0.25 

-0.32 

0.52 

0.39 

0.69 

0.64 

Near-bed  salinity 

0.07 

-0.52 

0.49 

0.35 

0.80 

0.72 

Surface  salinity 

0.11 

-0.61 

0.36 

0.02 

0.47 

0.24 

<j> 

0.07 

-0.62 

0.20 

-0.33 

0.36 

-0.03 

Table  4-1:  Squared  correlation  coefficients  (r  )  and  skill  scores  (SS)  for  the  North  flats.  (*)  One 
location  had  a  large  mean  difference  in  depth.  The  next  lowest  depth  skill  score  was  0.81 . 


The  model  simulations  are  compared  with  time  series  of  measured  water  depths,  cross- 
and  along-shore  velocities,  near-bed  and  surface  salinity  (which  dominates  the  water 
density),  and  stratification  (defined  by  the  potential  energy  anomaly,  O)  at  5  locations  on 
the  south-central  flats  in  June,  2009  (described  further  by  Ralston  et  al.,  2012),  and  at  25 
locations  on  the  north  flats  in  July,  2009  (red  and  yellow  symbols,  respectively,  in  Fig.  4- 
1).  Skill  scores  were  calculated  as  (Murphy  1988): 

S5  =  H-(r-J)2-(^)2  (1) 

where  r  is  the  correlation  coefficient,  and  sm  and  s(/  are  the  sample  standard  deviations 
and  m  and  d  are  the  sample  means  of  the  model  and  data,  respectively.  The  model 
simulates  accurately  the  depth  variations  (primarily  tidal  fluctuations)  (mean  r  and  SS  > 
0.9)  and  the  cross-shore  velocities  (mean  r  and  SS  >  0.5  on  the  north  flats  and  >  0.9  on 
the  south  flats).  The  alongshore  velocity  typically  is  weaker  than  the  cross-shore  velocity, 
and  is  more  sensitive  to  wind  and  pressure  gradients,  and  thus  more  difficult  to  simulate. 
Similarly,  surface  salinity  is  dependent  on  small-scale  structure  of  fronts  and  winds,  and 
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thus  it  is  not  simulated  as  well  as  near-bottom  salinity.  FVCOM  has  higher  skill  on  the 
south-central  flats  than  on  the  north  flats. 


South-central  flats,  June  2009 

Min 

Mean 

max 

r2 

SS 

r2 

SS 

r2 

SS 

Depth 

0.99 

0.99 

1.00 

1.00 

1.00 

1.00 

Cross-shore  Velocity 

0.88 

0.87 

0.93 

0.92 

0.95 

0.95 

Alongshore  Velocity 

0.39 

0.33 

0.50 

0.45 

0.62 

0.56 

Near-bed  salinity 

0.61 

0.42 

0.74 

0.59 

0.81 

0.72 

Surface  salinity 

0.40 

-0.01 

0.45 

0.17 

0.50 

0.36 

o 

0.40 

0.04 

0.44 

0.13 

0.52 

0.23 

Table  4-2:  Squared  correlation  coefficients  ( r  )  and  skill  scores  (SS)  for  the  South-central  flats 


The  potential  energy  anomaly  is  sensitive  to  both  surface  and  bottom  salinity,  and  thus 
predictions  are  sensitive  to  small  errors.  Compared  with  the  data  (Tables  4-1  and  4-2), 
FVCOM  underestimates  O,  which  may  be  partly  an  effect  of  the  reduction  of  maximum 
horizontal  density  gradients  by  finite  cell  sizes.  Although  the  correlation  coefficients 
between  observed  and  predicted  O  are  low,  the  FVCOM  model  exhibits  behavior 
consistent  with  the  data,  showing  freshwater  surface  plumes  spreading  over  the  flats  at 
variable  rates. 

4.1.4  Time  series  of  individual  terms  (phase-averaged). 

Model  results  for  the  phase-averaged  terms  of  the  O-balance  were  estimated  using 
simulated  water  depths,  velocities,  and  densities  at  the  locations  of  the  field 
measurements  discussed  in  Chapter  3.  Although  the  processing  for  model  and  data  was 
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similar,  the  available  model  results  are  for  July,  whereas  the  field  data  was  collected  in 


Figure  4-2:  Time  series  of  the  different  terms  of  the  phase-averaged  <I>  dynamic  equation  based 
on  measurements  (a  and  b)  and  FVCOM  model  results  (c  and  d)  for  type  1  (a  and  c)  and  type  2  (b 
and  d)  tides.  Measurements  were  collected  in  late  August  (as  in  Chapter  3)  and  model  results  are 
from  early  July.  The  shaded  region  indicates  the  water  depth.  FVCOM  results  are  for  the  location 
in  the  model  corresponding  to  the  measurement  location  in  Chapter  3. 

August.  Despite  the  different  time  periods,  the  stratification  balance  in  the  model  is 
qualitatively  similar  to  that  in  the  data  (Fig.  4-2).  Type  1  tides  (Fig.  4-3,  blue  shading) 
have  positive  depth-mean  straining  on  the  ebb  and  negative  depth-mean  straining  on  the 
flood  (Fig.  4-2a  and  c,  dashed  purple  curve).  On  type  2  tides  (Fig.  4-3,  red  shading)  data- 
and  model-based  depth-mean  straining  is  positive  during  most  of  the  prolonged  high 
water  (Fig.  4-2b  and  d,  dashed  purple  curve).  Mixing  (Fig.  4-2  dashed  green  curve)  is 
stronger  in  the  FVCOM  model  than  in  the  data,  but  tends  to  be  negative  on  ebb  tides  and 
during  the  prolonged  high  water  of  type  2  tides.  Advection  (Fig.  4-2,  solid  dark  blue 
curve)  is  smaller  in  the  FVCOM  model  than  the  data.  Similar  to  the  data,  modeled  non- 
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mean  straining  (Fig.  4-2  solid  light  blue  curve)  is  negligible  during  type  1  tides  but 
becomes  similar  in  size  to  the  other  terms  during  the  prolonged  high  water  of  type  two 
tides,  although  of  opposite  sign  to  the  data.  The  model  results  confirm  that  vertical 
advection  (Fig.  4-2d,  pink  dotted  curve)  partially  balances  the  non-mean  straining.  Depth- 
change  is  a  relatively  small  term  (Fig.  4-2,  dotted  yellow  curve). 


Yearday  2009 


Figure  4-3:  Simulated  water  depths  at  the  offshore  edge  of  the  north  flats  with  the  red  and  blue 
shaded  regions  indicating  type  2  and  type  1  tides,  respectively. 


Although  the  data-  and  model-based  processes  had  similar  magnitudes,  the  phase 
relationships  (e.g.,  the  timing  of  maxima  and  minima)  between  the  processes  are  different 
for  the  data  and  model,  presumably  owing  to  small  errors  in  location  and  frontal 
positions.  Consequently,  the  model  is  primarily  used  to  examine  the  magnitude  of  the 
tidal  variability  of  the  processes,  defined  as  the  rms  value  of  the  phase  averaged  process 
over  a  type  1  or  type  2  tidal  cycle,  rather  than  the  timing  of  the  intratidal  behavior. 

4.2  Results 

4.2.1  Process  activity  in  different  directions 

The  observations  (Chapter  3)  suggest  that  alongshore  processes  can  be  important  to  the 
stratification,  similar  to  ROFIs  and  other  regions  with  strong  alongshore  density  gradients 
and  flows.  Here  the  simulated  spatial  dependence  of  the  along-  and  across-flow 
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advection,  depth-mean  straining,  and  non-mean  straining  are  examined  during  type  1  and 
type  2  tides  (Fig.  4-3).  The  major  and  minor  axes  of  depth-mean  velocity  (Figs.  4-4c  and 
d),  rather  than  the  directions  along  and  across  bathymetric  contours,  are  chosen  to  define 
the  x-  and  y-coordinate  directions,  respectively,  because  it  is  the  mean  or  shear  velocity 
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Figure  4-4:  (a  and  b)  Ellipses  encompassing  the  major-  and  minor  axes  of  velocity  shear  over  a 
tidal  cycle  and  contours  of  tidally  averaged  depth-mean  salinity,  and  (c  and  d)  ellipses 
encompassing  the  major  and  minor  axes  of  depth-mean  velocity  over  a  tidal  cycle  and  contours  of 
tidally  averaged  O  for  type  1  (a  and  c)  and  type  2  tides  (b  and  d).  The  black  4x’  indicates  the 
central  location  examined  in  Chapter  3. 

acting  on  the  horizontal  density  or  <E>  gradient  that  produces  the  directional  stratifying 
processes.  Depth-mean  flows  tend  to  be  rectilinear  (and  aligned  with  the  cross-shore 
direction,  not  shown)  on  the  south  flats,  but  are  near  equal  in  both  directions  on  the  north 
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flats.  In  some  cases  the  velocity  major  axis  (x)  direction  is  roughly  alongshore  on  the 
north  flats. 

The  largest  tidal  variability  in  advection  and  depth-mean  straining  occurs  in  the  x- 
direction  in  the  gutter,  near  the  mouth  of  the  north  fork,  and  on  the  south-central  flats 
(dark  red  in  Figs.  4-5-  and  4-6a,  b,  d,  and  e).  The  maximum  rms  depth-mean  straining  on 
the  flats  occurs  slightly  closer  to  the  marshes  during  type  2  than  during  type  1  tides 
(compare  Fig.  4-6d  with  Fig.  4-5d).  Variations  in  advection  and  depth-mean  straining  on 
the  north  flats  are  about  a  factor  4  smaller  than  those  on  the  south-central  flats.  Non-mean 
straining  variability  typically  is  small,  except  in  the  gutter  (Figs.  4-5-  and  4-6g  and  h). 


The  advection  and  depth-mean  straining  have  more  variability  in  the  x-direction  than  in 
the  y-direction  at  most  locations  on  the  central  and  south  flats  (Figs.  4-5-  and  4-6c  and  f, 
greens  and  blues).  However,  although  the  maximum  y-direction  contribution  (-0.0006)  is 
about  a  factor  of  3  smaller  than  the  maximum  x-direction  variations  (0.0020),  the  rms 
advection  and  depth-mean  straining  along  the  y-axis  are  equal  to  or  greater  than  those  in 
the  x-direction  at  many  locations  on  the  north  flats  (Figs.  4-5-  and  4-6c  and  f,  pinks  and 
reds).  The  y-component  of  rms  depth-mean  straining  also  dominates  between  the 
channels  near  the  marshes  on  the  south-central  flats. 

The  importance  of  y-axis  terms  on  the  north  flats  is  related  to  the  large  tidal  variability  of 
the  depth-averaged  velocity  and  of  the  integrated  shear  velocity  along  the  minor  axis.  The 
latter  of  the  two,  integrated  shear  (Fig.  4-3a  and  b),  is  defined  as 


-fv  (u  -  u)-dz 

D J-hK  J  D 


(2) 


where  u  is  the  vector  velocity,  the  overbar  is  the  depth  average,  r\  is  the  surface  height,  h 
is  the  depth  of  the  bed,  D=r]+h  is  the  total  water  depth,  and  z  is  the  vertical  coordinate 
with  positive  upward.  Integrated  shear  is  used  rather  than  the  maximum  velocity 
difference  because  it  appears  in  the  DS  term  of  the  dynamic  O  equation  (Chapter  3, 
equation  (8)).  Shear  given  as  maximum  differences  in  velocity  is  between  10  and  100 
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times  as  large  as  the  integrated  shear,  and  typically  is  similar  in  magnitude  to  the  depth- 
averaged  velocity. 
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Figure  4-5:  Advection  (a,  b,  and  c),  depth-mean  straining  (d,  e,  and  f)  and  non-mean  straining  (g, 
h,  and  i)  for  type  1  tides.  The  columns  from  left  to  right  are  rms  value  in  the  x-direction  (a,  d,  g), 
rms  value  in  the  y-  direction  (b,  e,  h),  and  the  log  base  10  of  the  ratio  of  rms  y-axis  component 
divided  by  the  rms  x-axis  component  (c,  f,  i).  The  x-  and  y-directions  are  defined  by  the  velocity 
major  and  minor  axes  shown  in  Fig.  4-4c.  The  white  4x’  indicates  the  equivalent  of  the  central 
location  examined  in  Chapter  3. 
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Figure  4-6:  The  same  as  Fig.  4-5  but  for  type  2  tides.  The  x-  and  y-  directions  are  defined  by  the 
velocity  major  and  minor  axes  shown  in  Fig.  4-4d. 


On  the  south  and  central  flats,  the  tidal  velocity  is  nearly  rectilinear  (Fig.  4-4)  and  the 
major  axis  of  the  velocity  is  aligned  roughly  with  the  direction  of  the  density  gradient,  so 
rms  advection  is  dominated  by  the  x-component.  The  y-axis  component  of  straining  is 
slightly  larger  than  that  for  advection  on  the  south-central  flats  because  the  shear  is  not 
aligned  perfectly  with  the  mean  velocity  and  the  shear  ellipse  is  more  eccentric  than  that 
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for  mean  velocity  (compare  Figs.  4-4a  and  b  with  4-4c  and  d).  On  the  north  flats,  the 
density  and  O  gradients  are  smaller  and  not  as  well  aligned  with  velocity  or  shear  as 
those  on  the  south-central  flats  (Fig.  4-4),  and  thus  the  total  (x-  plus  y-direction)  rms 
variability  of  the  advection  and  straining  are  smaller.  However,  the  major  and  minor  axes 
of  the  velocity  are  nearly  equal  (Fig.  4-4),  and  the  relatively  strong  flows  and  shear  in  the 
minor  axis  direction  can  advect  and  strain  the  density  field  along  the  north  flats  leading  to 
relatively  large  ratios  of  y-  to  x-axis  contributions  (Figs.  4-5-  and  4-6c,  f,  and  i). 


Overall,  the  model  results  on  the  north  flats  show  increased  y-axis  contributions 
compared  with  the  field  data  results  shown  in  Chapter  3.  These  differences  are  partly  an 
effect  of  x-  and  y-  directions  being  defined  as  the  velocity  major  and  minor  axes  rather 
than  across  and  along  bathymetric  contours,  and  partly  a  result  of  different  velocity  and 
density  distributions.  Owing  to  this  difference,  it  is  necessary  to  interpret  the  model 
results  in  the  context  of  relative  distribution  of  velocity  and  density,  rather  than  directly 
applying  the  model  results  to  the  real  Skagit  Bay  tidal  flats.  Thus,  in  regions  with 
rectilinear  velocities  that  are  roughly  perpendicular  to  the  isopycnals  and  stratification 
contours,  such  as  near  the  south  fork  and  in  narrow  channels  or  estuaries,  the  processes 
creating  and  destroying  stratification  are  roughly  two-dimensional.  However,  in  regions 
with  a  relatively  large  velocity  contribution  in  the  minor  axis  direction  or  with  the  minor 
flow  axis  perpendicular  to  the  isopycnals,  such  as  can  occur  on  the  north  flats,  in  ROFIs 
(e.g.  de  Boer  et  al.,  2008),  and  on  estuarine  shoals  (e.g.  Ralston  et  ah,  2010a),  it  is 
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necessary  to  include  both  x-  and  y-components  to  estimate  the  effects  of  the  different 
processes  accurately. 
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Figure  4-8:  Cross-shore  distance  (positive  onshore)  versus  time  during  the  24-h  tidal  cycle  with 
contours  of  depth-averaged  total  salinity  gradient  magnitudes  for  phase-averaged  type  1  (a)  and 
type  2  (b)  tides.  Distance  is  along  the  transect  shown  in  Fig.  4-1. 


4.2.2  Differences  in  ebb  stratification  processes  between  type  1  and  type  2  tides 

The  maximum  rms  depth-mean  straining  on  the  flats  occurs  slightly  closer  to  the  marshes 
during  type  2  than  during  type  1  tides  (compare  Fig.  4-6f  with  Fig.  4-5f).  The  cross-shore 
transect  shown  in  Fig.  4-1  is  used  to  examine  this  difference  in  more  detail.  The  transect 
was  chosen  to  be  in  the  region  with  nearly  rectilinear  velocity  and  the  largest  straining 
signal  to  minimize  the  effects  of  y-direction  processes  and  to  focus  on  straining  effects. 
The  transect  profile  is  shown  in  Fig.  4-7.  The  more  shoreward  straining  during  type  2 
tides  is  owing  to  small  offshore  velocities  during  the  weak  ebb  allowing  the  salt  wedge, 
and  thus  the  region  with  the  strongest  horizontal  salinity  gradients  (Fig.  4-8),  to  remain 
close  to  the  marshes  for  a  longer  period  of  time.  In  addition,  during  type  2  tides, 
stratification  increases  throughout  the  weak  ebb  and  flood  owing  to  the  weak  flows  (and 
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mixing).  The  maximum  stratification  often  is  largest  following  the  weak  flood  (compare 
Fig.  4-9b  time=19  h  with  9a  time=7  h,  and  see  Ralston  et  al.,  2010a;  Giddings  et  al., 
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Figure  4-9:  Cross-shore  position  (positive  onshore)  versus  time  during  the  ebb  with  contours  of 
<I»  (a  and  b),  depth-mean  straining  (c  and  d)  and  mixing  (e  and  f)  for  phase-averaged  type  1  (a,  c, 
and  e)  and  type  2  (b,  d  and  f)  tides.  Times  are  subsets  of  those  in  Fig.  4-8. 


During  type  1  ebbs,  strong  velocity  shear  (not  shown)  and  moderate  horizontal  salinity 
gradients  (Fig.  4-8a)  cause  strong  straining  over  most  of  the  flats  throughout  most  of  the 
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O  (J/m3)  dO/df  (mW/m3) 


ebb  (yellow  to  red  contours  in  Fig.  4-9c),  while  mixing  remains  moderate  to  weak  owing 
to  moderate  velocities.  Thus,  the  stratification  increases  until  about  mid  ebb  (Fig.  4-9a, 
time  about  9  hrs).  During  type  2  ebbs,  O  is  large  and  increasing  until  time  ~  21  h  but  then 
weakens  rapidly  because  the  large  ebb  velocity  leads  to  mixing  outcompeting  straining 
(Figs.  4-9b,  d,  and  f).  Thus,  in  contrast  to  type  1  tides  in  which  the  water  on  the  offshore 
flats  and  the  gutter  is  stratified  at  low  tide  (Fig.  4-9a,  d>  ~  20  at  x  =  0  and  time  =12  h)  and 
de-stratifies  during  the  flood,  during  type  2  tides  the  water  on  the  flats  and  gutter  de- 
stratifies  on  the  ebb  (Fig.  4-9b,  d>  ~  0  at  x  =  0  and  time  =24  h)  and  the  flood  starts  with 
well-mixed  water. 
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Figure  4-10:  dO/dt  vs.  flat-averaged  instantaneous  Simpson  number  during  ebb  tides. 
Each  point  is  a  model  timestep  during  an  ebb  tide. 
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4.3  Discussion 


4.3.1  Ebb  tide  changes  in  stratification 


The  difference  in  stratification  magnitudes  at  the  end  of  the  type  1  and  2  ebbs  at  least 
partly  results  from  differences  in  the  ratio  of  straining  to  mixing.  This  relationship  can  be 
characterized  using  the  Simpson  number  (Si),  which  has  been  formerly  called  the 
horizontal  Richardson  number  (Burchard  et  al.,  2011;  Verspecht  et  al.,  2009;  Stacey  et 
al.,  2001;  Simpson  et  al.,  1990).  The  Simpson  number  can  be  conceptualized  by 
comparing  the  scaling  of  the  mixing  and  depth-mean  straining  terms  (Simpson  et  al., 
1990).  The  scaling  of  the  depth-mean  straining  term,  expressed  in  terms  of  gravity, 
horizontal  density  gradient,  velocity,  and  depth  is: 

[DS]*gdfxUD  (3) 

while  the  mixing  term  scaling  is: 

(4) 


Choosing  the  friction  velocity  as  the  appropriate  velocity  scale  after  the  scaling  shown  by 
Stacey  et  al.  (2001),  Si  can  be  expressed  as  the  ratio  of  [DS]  to  [M]: 


Si=^-2 

p o  dx  U£ 


(5) 


where  17*  is  the  friction  velocity,  D  is  the  water  depth,  x  is  horizontal  distance,  p  is  water 
density,  and  po  is  a  reference  density  (here  1000  kg/m  ).  For  the  following  calculations, 
all  quantities  are  averaged  across  the  tidal  flats  to  minimize  the  effects  of  advection  of  <t> 
and  it  is  assumed  that  on  average  the  density  gradient  and  the  depth-averaged  velocity  are 
aligned,  so  the  magnitudes  of  each  are  used  in  the  calculations.  The  resulting  tidal-flat- 
averaged  but  time-varying  value  of  Si  can  be  described  as: 


Si=JL|Vp|]£!i 
Po  rV.  I2 


where  the  vertical  bars  are  defined  as: 


(6) 


ItI 


£  XjAj 
X  A-i 
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where  x  denotes  an  arbitrary  spatially-varying  quantity,  Xi  is  the  value  of  that  quantity  at 
cell  or  node  index  i,  A;  is  the  area  associated  with  that  cell  or  node,  and  the  summation  is 
over  all  locations  on  the  tidal  flat  portion  of  the  domain  that  are  underwater. 


Figure  4-1 1 :  Local  d<b/dt  vs.  instantaneous  Simpson  number  during  ebb  tides.  The  location 
corresponds  to  that  examined  in  Chapter  3.  Each  point  is  a  model  timestep  during  an  ebb  tide  with 
a  high  (blue  circles)  or  low  (green  squares)  advection  term. 


For  Si>l,  d|<h|/dt  tends  to  be  positive,  indicating  increasing  average  O,  while  Si<l  is 
associated  with  negative  d|d>|/3t  (Fig.  4-10).  For  the  type  2  tides  (Fig.  4-10  green 
squares),  Si  is  less  than  unity  throughout  the  latter  half  of  the  ebb  and  the  water  becomes 
well  mixed.  For  type  1  tides  (Fig.  4-10  open  circles),  Si  is  rarely  less  than  unity,  and  the 
water  remains  stratified  at  the  end  of  the  ebb.  Decreasing  |d>|  does  occur  for  Si>l,  which 
is  a  result  of  the  mean  velocity  advecting  O  offshore  on  the  ebb,  off  of  the  tidal  flats 
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proper.  In  addition,  the  very  low  magnitudes  of  d|<h|/dt  occurring  at  small  Si  values  are 
associated  with  very  shallow  flows  where  the  water  column  is  already  mixed. 

Local  values  of  Vp,  I/*,  and  D  were  used  to  determine  a  local  value  of  Si  on  the  northern 
midflats  that  could  be  compared  with  the  data  presented  in  Chapter  3  (Fig.  4-11).  The 
assumption  of  parallel  density  gradient  and  flow  velocity  is  retained,  although  this  may 
overestimate  the  relevant  value  of  the  density  gradient  and  Si  in  the  north-western  region 
of  the  flats  (Fig.  4-4).  The  local  and  flat-averaged  d<b/3t  and  Si  (SC)  are  similar,  but  the 
local  values  are  more  often  in  the  lower-right  quadrant  of  high  SiL  and  negative  30/dt 
owing  to  the  effects  of  local  advection  overwhelming  the  balance  between  mixing  and 
straining  (Fig.  4-11).  Eliminating  data  points  corresponding  to  advection  values  higher 
than  the  median  removed  most  of  the  points  in  the  lower  right  quadrant.  The  remainder  of 
the  samples  in  the  lower  right  quadrant  may  violate  the  assumption  of  cross-isopycnal 
velocity.  Removing  times  with  high  advection  preferentially  removes  high  SiL  samples 
which  is  consistent  with  the  results  from  the  data  section  (Chapter  3)  showing  depth- 
mean  straining  and  advection  having  similar  magnitudes. 

Measurements  were  insufficient  to  determine  17*  directly  from  data,  so  a  proxy  based  on 
the  mean  velocity  was  used  to  calculate  SC: 

Ul  =  CdU 2  (8) 

The  value  of  Cd  used  for  the  data  was  0.0023,  which  was  determined  from  the 
relationship  between  I/*  and  U  in  the  FVCOM  model  at  the  midflat  location.  The 
FVCOM  model  correlation  between  I/*  and  U  was  similar  to  that  between  I/*  and  near¬ 
bottom  velocity.  The  field  data  show  that  higher  values  of  30/dt  were  associated  with 
higher  values  of  SC  (Fig.  4-12),  but  there  is  a  substantial  amount  of  scatter  that  is  not 
reduced  by  removing  points  with  high  values  of  advection  (not  shown).  The  scatter  may 
be  a  result  of  non-parallel  density  gradients  and  flows  at  the  measurement  location,  which 
is  neglected  in  this  formulation  of  SiL.  However,  a  larger  proportion  of  the  samples  from 
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the  data  than  the  model  had  high  SiL,  which  helps  explain  why  stratification  remained 
higher  in  the  data  than  the  model. 
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Figure  4-12:  Local  d<b/dt  vs.  instantaneous  Simpson  number  during  ebb  tides  using  field  data 
from  Chapter  3. 


Si  (and  variations  of  this  non-dimensional  number)  have  been  used  to  characterize  spatial 
and  temporal  variations  of  estuarine  dynamics  (e.g.,  Simpson  et  al.,  1990).  Here  it  is 
shown  that  an  instantaneous  Si  can  be  used  to  describe  transitions  from  increasing  to 
decreasing  O  on  ebb  tides.  However,  spatial  averaging  is  preferable  to  remove  advection 
effects  that  obscure  the  effects  of  mixing  and  depth-mean  straining.  Notably,  Si  is  useful 
in  tidal  flat  regions  like  this  one  where  the  tidal  depth  changes  are  comparable  to  the 
mean  water  depth,  and  not  only  in  regions  where  tidal  amplitudes  are  small  compared  to 
water  depth  (e.g.  Burchard  et  al.,  2011) 
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4.3.2  Density  effects  on  the  flow  structure 


The  FVCOM  model  was  also  used  to  investigate  the  effects  of  density  on  the  circulation. 
The  momentum  balance  in  FVCOM  uses  the  primitive  equations,  including  the 
hydrostatic  assumption,  and  can  be  written  as  (Chen  et  al.,  2003): 


du  ,  du  ,  du  ,  du  c 

—  +  U—  +  v—  +  w—  -  fv  = 

dt  dx  dy  dz 

dv  ,  dv  ,  dv  ,  dv  ,  r 

—  +  u  — +  v— +  w  —  +  fu  = 

dt  dx  dy  dz 


+£  (*-»£) +F« 


dP_ 

dz 


=  -99 


(9) 

(10) 

(11) 


where  u,  v,  and  w  are  the  velocity  in  the  x,  y  and  z  directions,  and  z  is  vertical  up,  t  is 


time, /is  the  Coriolis  parameter,  p  is  density  and  po  is  a  reference  density,  P  is  pressure, 


Km  is  the  eddy  viscosity,  Fu  and  Fv  are  horizontal  momentum  diffusion  terms,  and  g  is 
gravitational  acceleration.  Non-uniform  density  affects  the  first  two  terms  on  the  right 
hand  side:  the  pressure  gradient  term  and  the  shear  or  eddy  viscosity  term,  respectively. 


65 


E 


O) 

c 


n 

o 


5356 

5352 

5348 

5344 


530  534  538  542  546 


Easting  (km) 


Figure  4-13:  Ellipses  encompassing  the  major  and  minor  axes  of  depth-averaged  velocity  for  the 
base  (blue)  and  no-density  cases  (red).  The  pink  circle  represents  0.2  m/s. 


Integrating  the  vertical  pressure  gradient  over  the  water  column  results  in  pressure  as  a 
function  of  depth: 

P(z)  =  P0  + £  pgdz  (12) 

where  Po  is  atmospheric  pressure  (uniform),  and  r/  denotes  the  level  of  the  surface. 
Taking  a  horizontal  derivative  of  P(z)  and  using  the  Leibniz  integral  rule  results  in: 
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(13) 


dP(z) 

dx 


dr\  rr\  dp  , 

=  P,9y,  +  Sz  T,3dz 


Z  dx  ‘ 


where  pv  is  the  water  density  at  the  surface.  The  first  term  on  the  right  is  the  barotropic 
pressure  gradient,  and  the  second  term  is  the  baroclinic  pressure  gradient,  which  is 
dependent  primarily  on  horizontal  density  distributions,  but  also  may  be  modified  by  the 
vertical  density  distribution  (or  stratification).  The  baroclinic  pressure  gradient  will  drive 
flows  from  regions  of  heavier  to  lighter  water.  The  force  associated  with  the  baroclinic 
pressure  gradient  will  tend  to  be  small  near  the  surface  and  increase  in  magnitude  to  the 
bed,  but  will  be  all  in  the  same  direction  given  a  uniform  density  gradient  and  can  thus 
affect  both  the  depth-averaged  flows  and  the  shear. 


Stable  stratification  leads  to  conversion  of  turbulent  energy  to  potential  energy  as  water 
parcels  move  vertically  (rather  than  transfers  of  momentum),  which  is  expressed  as  a 
reduction  in  the  eddy  viscosity.  It  can  also  inhibit  the  generation  of  turbulence  in  the  first 
place.  Overall,  stratification  results  in  greater  total  shear  in  the  water  column  (Turner, 
1973)  as  there  is  less  interaction  between  water  layers  and  the  bottom  stress  slows  the 
nearbed  water  more  than  the  surface  water.  The  effect  of  stratification  on  eddy  viscosity 
is  similar  to  that  on  eddy  diffusivity  (Kv)  for  which  one  parameterization  is  (Munk  and 
Anderson  1948): 

Kv  =  K0(  1  +  3.33i?i)-i  (14) 

where  Ko  is  eddy  diffusivity  without  stratification  and  R,  is  the  Richardson  number, 
which  increases  with  increasing  vertical  density  gradients.  Stratification  will  thus  always 
affect  the  distribution  of  the  velocity  in  the  water  column.  Stratification  in  the  absence  of 
horizontal  density  gradients  may  also  affect  the  depth-averaged  velocity  through  a 
reduction  in  the  bottom  drag  coefficient,  requiring  greater  water  velocity  to  acheive 
bottom  stress  that  balances  the  barotropic  pressure  gradient. 
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Figure  4-14:  Residual  flow  vectors  for  the  base  (blue)  and  no-density  cases  (red).  The 
pink  arrow  corresponds  to  0.1  m/s. 

To  illustrate  the  effects  of  density  and  stratification  on  the  flows,  The  FVCOM  model  is 
run  with  the  same  initial  conditions  (including  spatial  salinity  gradients),  but  with  the 
salinity  treated  as  a  passive  tracer,  removing  all  effects  of  density  differences  (no-density 
case).  Thus,  in  the  no-density  case,  the  baroclinic  pressure  term  and  the  density  effect  on 
the  eddy  viscosity  are  neglected.  These  results  are  compared  with  a  standard  (base  case) 
model  run  for  late  June.  Major  and  minor  axes  of  velocity  and  integrated  shear  are 
determined  for  each  two-week-long  model  run. 
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Figure  4-15:  Ellipses  encompassing  the  major  and  minor  axes  of  depth- integrated  shear 
for  the  base  (blue)  and  no-density  cases  (red).  The  pink  circle  represents  0.02  m/s. 


The  variability  over  a  tidal  cycle  of  the  depth-averaged  velocity  is  similar  for  both  model 
runs  (Fig.  4-13).  The  largest  differences  in  the  ellipses  shown  for  mean  velocity  occur  on 
the  offshore  flats  and  in  the  gutter.  The  baroclinic  pressure  gradient  scales  with  the  depth 
as  well  as  the  density  gradient,  and  thus  tends  to  be  smaller  in  shallow  water.  However, 
although  the  depth-averaged  baroclinic  pressure  gradient  is  small,  there  is  a  measureable 
effect  on  the  magnitude  and  direction  of  the  residual  flows  (Fig.  4-14).  Furthermore,  the 
density  gradients  and  stratification  significantly  increase  the  shear  (Fig.  4-15),  especially 
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in  the  gutter  and  on  the  north  flats  where  the  tides  are  not  rectilinear.  Even  in  the  shallow 
areas  by  the  south  fork,  the  variable  density  makes  a  difference  in  shear.  Thus, 
stratification  will  more  significantly  affect  processes  that  depend  on  the  depth  distribution 
of  the  velocity,  such  as  those  driven  by  nearbed  flows. 

4.4  Summary 

On  the  south-central  flats,  advection  and  straining  are  roughly  aligned  with  the  velocity 
major  axis,  which  is  roughly  aligned  with  the  density  gradient  direction.  In  contrast,  on 
the  north  flats,  the  major  and  minor  axes  of  the  tidal  velocity  variability  are  nearly  equal 
but  neither  is  aligned  with  the  density  or  O  gradient,  and  the  advection  and  straining  have 
nearly  equal  contributions  along  the  major  and  minor  axes  of  the  velocity.  Likewise, 
observations  on  the  north  flats  (Chapter  3)  showed  that  velocity  magnitudes  were  similar 
in  all  directions,  so  the  orientation  of  the  processes  depended  on  the  density  gradient, 
which  was  primarily  cross-shore.  In  addition,  the  FYCOM  model  shows  stratifying 
processes  that  are  stronger  overall  on  the  southern  flats  because  of  the  near-alignment  of 
velocity  or  shear  and  density  gradient  in  that  location. 

The  importance  of  mixing  to  destratifying  the  water  column  depends  on  the  tidal  range. 

In  this  early  July  period  type  1  tides  had  a  relatively  small  tidal  range,  and  stratification 
persists  until  low  tide  because  the  velocities  remain  relatively  weak  and  mixing  is  small. 
In  contrast,  the  early  July  type  2  tides  with  a  large  total  tidal  range,  and  strong  flows  mix 
the  water  column  during  the  strong  ebb.  The  Simpson  number,  which  takes  into  account 
the  horizontal  density  gradient,  the  depth,  and  the  friction  velocity,  can  be  used  to 
evaluate  the  relative  importance  of  mixing  and  straining.  Mixing  dominates  over  straining 
when  the  instantaneous  flat-averaged  or  local  Simpson  number  is  less  than  unity,  and 
straining  dominates  when  the  Simpson  number  is  greater  than  unity.  Local  advection 
causes  scatter  in  the  local  results.  The  Simpson  number  also  was  useful  in  determining 
the  relative  importance  of  mixing  and  straining  to  the  field  data. 
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Although  the  tidal  variability  of  depth-averaged  velocity  is  insensitive  to  density  effects, 
the  variability  of  the  velocity  shear  and  the  magnitude  and  direction  of  the  residual  flows 
are  affected  by  the  inhomogeneous  density. 
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5  Conclusions  and  future  directions 


In  this  thesis,  observations  and  FYCOM-simulations  of  water  levels,  velocities,  and 
densities  in  Skagit  Bay,  Washington,  are  used  to  investigate  the  processes  leading  to 
production,  destruction,  and  movement  of  stratification  on  the  tidal  flats.  Conclusions  for 
each  of  the  primary  questions  are  listed  below. 

1 .  What  processes  affect  the  stratification,  and  what  causes  spatial  variations  in  the 
magnitudes  of  these  processes? 

The  dominant  processes  controlling  stratification  on  the  Skagit  Bay  tidal  flats  are 
advection  and  depth-mean  straining,  with  mixing  on  strong  ebb  tides  and  during  periods 
of  high  stratification.  On  the  strong  floods,  a  salt  wedge  with  a  near  vertical  interface  is 
transported  onshore,  while  on  strong  ebbs  the  salt-wedge  is  strained  and  mixed  as  it  is 
advected  offshore.  The  magnitude  of  the  advection  (depth-mean  straining)  depends  on  the 
magnitude  of  the  stratification  (density)  gradient  and  the  velocity  (shear)  as  well  as  the 
angle  between  them.  Numerical  simulations  (using  FYCOM)  suggest  that  stratification  is 
stronger  on  the  south  flats  because  the  velocity  major  axis  and  the  density  gradients  are 
roughly  aligned.  These  results  suggest  that  narrow  estuaries  and  tidal  flat  channels  where 
the  major  flow  axis  is  aligned  with  the  salinity  gradient  are  likely  to  have  larger 
magnitudes  of  stratification  (for  fixed  freshwater  discharge  and  tidal  velocities)  than 
regions  with  similar  flow  magnitudes  in  both  directions  and  strong  along-isopycnal  flows. 

2.  What  is  the  relative  importance  of  cross-  and  alongshore  processes,  and  what  controls 
this  balance? 

On  the  south  flats,  tidal  velocities  are  roughly  rectilinear,  and  thus  processes  are  small  in 
the  minor  axis  direction  (roughly  alongshore).  On  the  north  flats,  the  shear  is  roughly 
parallel  with  the  isopycnals,  and  thus  there  is  little  alongshore  straining  despite  strong 
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alongshore  shear  during  the  type  2  higher  low  water,  except  on  the  offshore  portions  of 
the  north  flats.  However,  alongshore  advection  is  significant  over  most  of  the  north  flats 
because  stratification  is  spatially  variable  and  alongshore  flows  are  strong  during  the 
weak  ebb  and  flood.  Owing  to  the  alongshore  advection,  cross-shore  straining,  and  weak 
depth-averaged  flows  (resulting  in  little  mixing)  during  this  period,  stratification 
increases  until  the  second  high  tide,  resulting  in  larger  magnitudes  of  stratification  during 
type  2  (nearly  diurnal)  than  type  1  (roughly  semi-diurnal)  tides,  similar  to  results  from 
prior  studies  of  salt-wedge  estuaries. 

3.  Where  is  stratification  generated  and  destroyed? 

The  observations  suggest  that  stratification  is  generated  near  the  measurement  location  on 
the  north  flats  (straining  dominates  over  mixing),  which  is  about  1.5  m  shallower  than 
lower  low  water.  Decreased  stratification  often  was  owing  to  advection,  as  the  stratified 
water  is  transported  towards  the  central  flats,  or  out  of  the  bay  and  into  the  deeper  water 
of  Puget  Sound  and  the  Strait  of  Juan  de  Fuca.  Numerical  (FVCOM)  simulations  suggest 
the  dominance  of  generation  (owing  to  straining)  or  destruction  (owing  to  mixing)  of 
stratification  on  the  flats  reflects  an  instantaneous  Simpson  number  greater  or  less  than 
unity,  respectively. 

Future  work  related  to  these  results  includes: 

•  Comparing  model  predictions  of  spatial  variability  of  stratification  for  low  and 
high  discharge  periods. 

•  Quantitatively  comparing  field  data  results  from  different  time  periods 

o  Although  the  combined  cross-  and  alongshore  transects  used  here  were 
deployed  only  during  late  August,  separate  cross-  and  alongshore  transects 
were  deployed  from  late  June  until  mid- August  (including  periods  with  a 
range  of  river  discharge  values). 

•  Evaluating  the  effects  of  wind  and  wave  forcing. 
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Exploring  the  relevance  of  other  non-dimensional  numbers. 


Future  work  that  would  improve  our  understanding  of  hydrodynamics  on  sandy  tidal  flats 
with  freshwater  input  includes: 

•  Collecting  field  data  to  evaluate  stratification  and  mixing  processes  during  storms 
with  significant  rainfall,  winds,  and/or  waves. 

•  Using  model  simulations  and  new  field  data  to  examine  the  effects  of  the  presence 
of  multidirectional  processes  on  the  sediment  transport  near  the  north  flats 
(including  trapping  of  sands  on  the  Skagit  flats). 

•  Comparing  Skagit  Bay  to  flats  on  the  sides  of  estuaries  to  examine  the  effects  of 
different  orientations  of  the  fresh-  and  saltwater  sources 


75 


76 


Appendix:  Mixing  Parameterization 

The  mixing  term  M  in  the  stratification  balance  uses  a  common  parameterization  of  the 
vertical  buoyancy  flux  (Nepf  and  Geyer,  1996;  Becker  et  al.,  2009;  Ralston  et  al.,  2010b) 
based  on  the  eddy  diffusivity  Kv.  The  eddy  diffusivity  is  estimated  using 
parameterizations  of  the  turbulence  produced  in  the  bottom  boundary  layer  modified  by 
the  ability  of  stratified  flow  to  support  turbulence  (Munk  and  Anderson,  1948).  The  eddy 

diffusivity  was  calculated  by  Kv  =  A',,)!  +  3.33RJ  ,  where  K0  is  the  estimated  eddy 

diffusivity  for  an  unstratified  water  column  and  R,  is  the  bulk  Richardson  number. 
Extending  a  boundary  layer  method  (Becker  et  al.,  2009),  K0  is  calculated  by 

/  z  \ 

K0  =  OAC0u*zb  1  -  —  where  Co  =  0. 1  is  based  on  a  fit  to  the  stratification  balance  and  to 

\  D) 

model  predictions  of  Kv  (Ralston  et  al,  2012)  u*  is  the  friction  velocity,  zb  is  the  smaller 
of  the  height  above  the  bed  or  the  thickness  of  the  bottom  boundary  layer,  and  D  is  the 
total  water  depth.  Co  was  added  here  as  a  scaling  factor  since  the  original  boundary  layer 
parameterizations  were  for  well-mixed  or  weakly  stratified  conditions  and  strongly 
overpredict  the  eddy  diffusivity  for  the  highly  stratified  conditions  found  on  the  Skagit 
Bay  tidal  flats.  An  isotropic  estimate  is  used  for  the  friction  velocity,  given  by  u2  = 
Cd(ii2  +  v2),  where  u  and  v  are  the  depth-averaged  cross-  and  alongshore  velocities, 
respectively,  and  the  drag  coefficient  Q  is  0.001.  The  data  collection  methods  did  not 
allow  for  more  accurate  estimates  of  the  friction  velocity  by  direct  measurement  of 

gD  A  p 

turbulence.  The  Richardson  number  is  calculated  by  Rt  =  ^ — ~ — -r ,  where  g  is 

p  (Aw)  +(Av) 

gravitational  acceleration,  p  is  the  depth-averaged  density,  and  Ap ,  Am,  Av  are  the 
maximum  difference  in  density,  cross-shore  velocity,  and  alongshore  velocity, 
respectively  (Byun  and  Wang,  2005).  The  thickness  of  the  bottom  boundary  layer  ( HBbl ) 

/  t-V>£ 

g  dp 

is  estimated  by  H BBL  =  m,  — - (Stacey  and  Ralston,  2005). 

\pRt  dx) 
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Alternate  methods  evaluated  for  the  mixing  parameterization  include  a  constant  and 

g  dpt  6uj  \2  <5v.  \2 

uniform  Kv,  a  depth- varying  gradient  R,  calculated  by  a,  =  —  — —  — —  +  —  ,  a 

p  oz  [\  oz  )  \oz  J 

constant  and  uniform  Kq,  a  depth-uniform  K()  calculated  by  K0  =  0.4 «  and 

combinations  thereof.  Unrealistic  values  of  Rj  are  obtained  when  calculated  in  a  depth- 


varying  manner  owing  to  the  constant  extrapolation  method  used  to  obtain  density  and 
velocity  values  throughout  the  water  column.  Constant  and  uniform  K0  or  Kv  were 
rejected  because  they  do  consider  the  effects  of  stratification  and  flow  speed.  The  choice 
of  the  method  used  was  based  on  the  best  fit  of  the  total  dynamic  balance. 
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